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Preface 



This book contains the written contributions to the International Symposium on 
Medical Simulation (ISMS ’04) held in Cambridge, Massachusetts, USA on June 17 th 
and June 18 th , 2004. 

Manuscripts are organized around five thematic sections relating to the 
multidisciplinary field of Medical Simulation: Soft Tissue Properties and Modeling, 
Haptic Rendering, Real-Time Deformable Models, Anatomical Modeling, and 
Development Frameworks. 

The objectives of the symposium are to gather researchers to present their most 
recent, and promising work, to highlight research trends and foster dialogue and 
debates among participants. Live demonstrations are also included at the meeting, but 
cannot be included in this volume. Finally, to address questions about areas for 
improvement and future directions of the field, we organized a panel of experts 
including technical, medical and educational representatives. 

This event continues the successful symposium organized by Herve Delingette and 
Nicholas Ayache, in France in June 2003. At that meeting we agreed that it would be 
beneficial for the community to have an annual gathering for the medical simulation 
community where researchers can exchange ideas and share their work in this 
emerging field. ISMS ’04 is co-organized by CIMIT / Harvard Medical School and 
Rutgers University. 

We received 50 submissions from 14 different countries. Each was evaluated by 
three members of the scientific committee. We selected 16 manuscripts for oral 
presentation and 16 for poster presentation. All accepted manuscripts were allowed a 
written contribution of equal length. These published contributions are from research 
institutes, universities or companies from our diverse research community, including: 
Germany, Russia, Netherlands, France, Switzerland, Belgium, Spain, United 
Kingdom, Italy, Japan, Korea, Israel, Canada and the United States of America. 

The quality of the selected contributions implies that the conference will be an 
important milestone in the development of this new, but rapidly growing field at the 
convergence of several disciplines with key applications for the future of health care. 
We welcome all participants to this intense and stimulating scientific event 



April 7, 2004 



Stephane Cotin and Dimitris Metaxas 




Introduction 



Medical Simulation - The Second Phase 



It has been well over a decade since medical simulation began moving from the 
laboratory to resident and student training. Numerous simulation companies have 
emerged and disappeared, academic studies have begun the validation process, and 
opponents now have to reconsider their standpoint. It appears this initial work in 
simulation is leading to a second phase - a phase that will take simulators out of the 
laboratory and into daily training. 

Many factors account for this transition. The technology has improved to a point 
where low cost systems can now achieve reasonably high fidelity. 

Initial validation studies demonstrate unequivocally that training on a simulator 
improves technical skills proficiency, performance in the operating room, reduces 
time and variability in performance, and helps eliminate errors. 

Training curricula are being developed to encompass all types of simulators, 
focusing the training on comprehensive curricula rather than simply upon the 
simulator. Current validation studies prove that high visual realism is not required to 
provide excellent training assessment and transfer. At times novices prefer to have a 
more abstract and simple simulation to correctly train in the more fundamental skills. 
Program directors, medical societies and certification boards have acknowledged the 
importance of simulation for training and assessing competence and are now actively 
supporting this new field. After the recent United States federal regulation requiring 
residents to limit work hours to 80 per week, participants in resident training are 
feeling the pinch. The potential of simulation to provide some of the training that is 
lost with this regulation is significant. 

Now is the time to incorporate simulation as an official part of both training and 
assessment of competency. Having made that statement, there are only a handful of 
simulators currently available. An enormous amount of hard labor is required to build 
new simulators, validate them and incorporate them into rigorous curricula for 
training, assessment and eventually certification. Medical simulation has finally 
reached the plateau which flight simulators reached in the 1950's, when the FAA 
officially required them for flight training and certification. Flight simulation 
required tens of billions of dollars and over five decades to achieve the extraordinary 
level of fidelity and reliability that it enjoys today - we can expect no less an 
investment in time, money and talent to reach the same level of excellence for 
medical simulation. 
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Experimental Observation and Modelling 
of Preconditioning in Soft Biological Tissues 



Alessandro Nava 1 , Edoardo Mazza 1 , Oliver Haefner 1 , and Michael Bajka 2 



1 Centre of Mechanics, ETH Zurich, 8092 Zurich, Switzerland 
nava@imes .mavt . ethz . ch 
Tel. (+41) 1 6327755 

2 Department of Gynecology, University Hospital Zurich, Switzerland 



Abstract. Constitutive models for soft biological tissues and in particu- 
lar for human organs are required for medical applications such as surgery 
simulation, surgery planning, diagnosis. In the literature the mechanical 
properties of biosolids are generally presented in “preconditioned” state, 
i.e. the stabilized conditions reached after several loading-unloading cy- 
cles. We hereby present experiments on soft tissues showing the evolution 
of the mechanical response in a series of loading and unloading cycles. 
The experimental procedure applied in this study is based on the so 
called “aspiration experiment” and is suitable for in-vivo applications 
under sterile conditions during open surgery. In the present study this 
technique is applied ex-vivo on bovine liver. A small tube is contacted 
to the target organ and a weak vacuum is generated inside the tube ac- 
cording to a predefined pressure history. Several identical loading and 
unloading cycles are applied in order to characterize the evolutive be- 
haviour of the tissue. The experimental data are used to inform the 
fitting of uniaxial and threedimensional continuum mechanics models. 
This analysis demonstrates that a quasi-linear viscoelastic model fails 
in describing the observed evolution from the “virgin” to the precondi- 
tioned state. Good agreement between simulation and measurement are 
obtained by introducing an internal variable changing according to an 
evolution equation. 



1 Introduction 

The mechanical characterization of soft biological tissues is a key issue for a 
large number of medical applications, such as surgery planning, surgical training 
deploying virtual reality based simulators, or diagnosis (see [1] , [2] ,[3] , [4] ) . Quan- 
titative sets are available on the mechanical properties of soft tissues, however 
very limited data are available on the in-vivo behaviour of soft tissues associated 
with human organs ([5], [6], [7], [8]). 

When subjected to repeated loading and unloading cycles the mechanical 
response of soft tissues changes significantly in early cycles and than stabilizes. 
This tendency to approach a stable condition in repeated cycles is called pre- 
conditioning, [9]. Fung [9] proposes to focus the attention to the mechanical 
response in the preconditioned state. This approach enables the use of the so 
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called “quasi linear viscoelastic” model. However, such a model might not be 
suitable for medical applications. In fact, the mechanical response of human or- 
gans during surgery correspond to a virgin state rather than a preconditioned 
state. 

The evolution of tissue behaviour in preconditioning is related to changes 
in its internal structure [9]. Investigations of the preconditioning behaviour can 
therefore provide relevant information on the capabilities of the tissue to adapt 
to a mechanical load and eventually to recover its original properties if unloaded. 
These characteristics might be useful in assessing the health status of the organ 
under consideration. 

We present here an experimental procedure for investigating the precondi- 
tioning behaviour of soft tissues. The results of ex- vivo cyclic tests with bovine 
liver are analyzed. The experimental data are rationalized by introducing a his- 
tory dependent state variable into the quasi-linear viscoelastic model. The inter- 
nal variable is determined by an evolution equation for its time rate of change. 

2 Experiments 

The aspiration test device shown in figure 1 has been developed by V. Vuskovic 
[5] . The device has been designed for in-vivo applications addressing issues asso- 
ciated with: safety, sterilizability, space limitation and a short data acquisition 
cycle time. Several modifications in the acquisition and in the control system 
have been introduced in order to make the device suitable for the current re- 
search requirements. 

The principle of working of the device is based upon the pipette aspiration 
technique [10]. The device consists of a tube in which the internal pressure can 
be controlled according to a desired pressure law. The investigation is performed 
by (i) gently pushing the tube against the tissue to ensure a good initial contact, 
(ii) creating a (time variable) vacuum inside the tube so that the tissue is sucked 
in through a smooth edged aspiration hole (diameter of 10 mm), see figure 3. 

Assuming the tissue to be isotropic and homogeneous in the small portion 
under deformation, a complete description of the deformed tissue can be given 
by simply monitoring the side-view profile of the tissue during its deformation. 
An optic fiber connected to an external source of light provides the necessary 
illumination in the inner part of the tube. 

The images of the side-view (figure 2) are reflected by a mirror and are 
captured by a digital camera mounted on the upper part of the device. The 
grabbed images are processed off-line in order to extract the profiles of the 
deformed tissue. 

The pressure inside the device is controlled by means of a pump, an air 
reservoir and two valves. With such a control system, a predefined pressure law 
can be accurately realized. In particular a pressure history of identical repeated 
cycles is imposed to the tissue in the present experiments: eight identical loading 
cycles are performed for a total duration of about 200 seconds, see figure 3. 

Time histories of measured pressure and deformation profiles are the input 
data used to evaluate the mechanical properties and to determine the constitu- 
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Fig. 1. Aspiration device and principle Fig. 2. Image grabbed by the digital 
of working. camera. 




Fig. 3. Pressure law imposed: 8 loading 
cycles. 
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Fig. 4. Displacement of the middle 
point of the “bubble” . 



five model. Fig. 4 shows the displacement of the middle point of the “bubble” 
as function of the time. 

Material softening is clearly visible from the increasing displacement mea- 
sured cycle after cycle. 

The experiments were performed ex- vivo on bovine liver, approximately 6 — 7 
hours after extraction. Since the aspiration technique does not require cutting 
of samples, the organ was tested as a whole organ. To prevent hydration, the 
organ was laterally immersed in a saline solution. No water drops were detected 
on the surface of the organ after the aspiration experiment was performed. 

3 Mechanical Modelling 

3.1 Uniaxial Model 

First evaluations of the experimental data obtained were performed by modelling 
the experiment as a 1-D problem, simply by correlating the pressure imposed 
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Fig. 5. 1-D models investigated. 




Fig. 6. Model a): Voigt body + Maxwell Fig. 7. Model b): Voigt body, composed 
body, composed of linear springs and of a linear dashpot and a non-linear 
daslrpot. spring (5th order polynomial law). 



on the surface of the tissue and the displacement of the highest point of the 
aspirated tissue “bubble”. This approach enables to identify the main physical 
features of the problem without the complications of a three dimensional model. 

Uniaxial models of linear and quasi-linear viscoelasticity (fig. 5), based on 
combination of linear dashpots and linear (a) and non linear (b) springs, were 
used. The model parameters were obtained by matching calculation and mea- 
surement for the first loading cycle. The comparison between calculated and 
measured curves are shown in fig. 6 and fig. 7. 

These results demonstrate the problems associated with quasi-linear vis- 
coelastic models: the use of a non-linear spring leads to a better description of 
the observed behaviour, but softening cannot be described with this approach. 
An extension of the uniaxial models of fig. 5 as a series of Voigt models has been 
evaluated and shown not to significantly improve model prediction. 

The following modification to the non linear Voigt model ((b) in fig. 5) is 
introduced in order to describe material softening: a non linear spring is used 
in which the parameters describing the current spring response ( K n ) are related 
to the corresponding constants in the virgin state /\ n0 by means of a softening 
variable (SV) which is function of the deformation history ( SV = 0 at t = 0) 



1 + SV 



(1) 



Experimental Observation and Modelling of Preconditioning 



5 




U, D2 

+1.486e+00 
■fl.363e4.00 
■f 1 . 239e-f00 
•fl.ll5e-f00 
+ 9 . 910e-01 
♦8 . 671e-01 

+7 . 432e-01 

♦ 6 . 194e-01 
4-4 . 955e-01 

+3 . 716e-01 

♦ 2 . 477e - 01 

♦ 1 . 239e-01 
-2 . 544e-22 





c 



B 



Fig. 8. Modified model b): Voigt body Fig. 9. Finite element model, 

with a non linear spring and the soften- 
ing variable SV. 



where the evolution equation for the time rate of change of SV is given by 

SV = a\u\ 0 (2) 

in which u is the displacement of the middle point of the “bubble” and a and (3 
are material constants. 

Fig. 8 shows the results with this new model: with the softening variable it is 
possible to obtain a good representation of the first cycle as well as the evolution 
over repeated cycles. 

The tissue is expected to return to the original state after some times without 
external loads. A recovery term such as the one used in the model by Rubin and 
Bodner [11] should be included in order to describe this effect. The present model 
does not describe recovery, since the experiments do not provide information for 
fitting a recovery term. These results provided the basic idea for the formulation 
of the three dimensional constitutive model. 



3.2 Three Dimensional Constitutive Model 

A three dimensional constitutive model is required for predicting the mechanical 
behaviour in structures (organs) of arbitrary geometry subjected to arbitrary 
kinematic and kinetic boundary. The quasi-linear viscoelastic model previously 
used for evaluating the aspiration experiment [12] has been modified according 
to the idea of an internal variable SV derived in the previous section. 

We consider the material as a homogeneous, isotropic continuum and due to 
the high water content, as incompressible. 

Hyperelastic materials are described in terms of a strain potential energy U 
which defines the free energy of the material as a a function of the deformation. 
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We used the so called reduced polynomial form [13]. 

N 

U = Y.C n (h-Z) n ( 3 ) 

n—1 



where 1 i = Af + A| + Ag is the first strain invariant (A;: principal stretches). 

Viscoelasticity is taken into account by applying time dependent (relaxation) 
coefficients to the constants that define the energy function: 

C n (t) = C n0 (l-Yj£(}-e-%) ) J (4) 



where C n o describe the instantaneous elastic response and Tj^ and Tk characterize 
the relaxation behaviour. In the current implementation we chose to stop at the 
fifth order of the series expansion for the strain potential energy (N = 5) and at 
the fourth order for the Prony series ( K = 4). 

The softening of the material is modeled by relating the coefficients C n o to 
the deformation history: 



C n o — 



a 



nO 



1 + SV 



( 5 ) 



where C n o are the elastic parameters of the “virgin” material and SV is deter- 
mined by the following evolution equation (SV = 0 at t = 0): 



SV = /J (h - 3) 



( 6 ) 



where /i is a material constant. 

The experiment is simulated by an axisymmetric finite element (FE) model, 
fig. 9. The FE program Abaqus 6.2 has been used for this purpose [14]. The 
starting geometry of the tissue is the one grabbed by the camera at the beginning 
of the experiment, at atmospheric pressure inside the tube. The force by which 
the tube is pressed against the soft tissue leads to a non-zero initial deformation. 

The FE calculation proceeds by applying the measured pressure history to 
the free tissue surface (D-X, Fig. 9). Axisymmetric hybrid triangles are used with 
linear pressure and quadratic displacement formulation. The overall dimensions 
of the model are selected in order to minimize the influence of the boundary 
conditions at the bottom (A-B) and at the side (B-C) on the displacement of 
the aspirated tissue “bubble”. The contact between the tissue and the device 
(C-D) is modelled as rigid-deformable contact with sliding. 

Similar as for the uniaxial models, the 13 material constants (C n o, and 

Tk) related to the lryperelaticity and viscoelasticity are determined from the first 
loading cycle. The constant related to the softening (/z) is determined in a second 
step by running a new optimization on the whole experiment (8 loading cycles). 

The optimization algorithm aims at minimizing the error function E: 



e = J2 

i 



I Zi - Zi I 



Zi 



( 7 ) 
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where z% and 2, are the measured and the calculated displacement of the point 
X respectively. The implemented optimization procedure is the Nelder-Mead 
simplex (direct search) method [15]. 

Fig. 10 and fig. 11 compare the predictive capabilities of a quasi-linear visco 
elastic model and a model which includes the softening variable SV respectively. 
The evolution of the mechanical response can be reproduced using the internal 
variable SV. 





experiment 

- - FEM with softening 



0 20 40 60 80 100 120 140 160 180 200 

Time [s] 



Fig. 10. Calculated (dashed line) and 
measured (continuous line) vertical dis- 
placements using the classic quasi-linear 
viscoelasticity. 



Fig. 11. Calculated (dashed line) and 
measured (continuous line) vertical dis- 
placements using the proposed 3-D 
model with the softening variable. 



4 Conclusions and Outlook 

An experimental technique, based on the aspiration experiment, for studying 
the behaviour of a tissue subjected to repeated cycles has been described. The 
control system of the aspiration device allows defining an arbitrary pressure 
history and performing a relatively large number of identical loading cycles. The 
procedure for cyclic experiments can be applied in-vivo during open surgery for 
testing the preconditioning behaviour of soft tissues in human organs. 

The analysis of the experimental results has shown that a quasi-linear vis- 
coelstic model is not adequate for describing the material behaviour in repeated 
cycles. A modification of the quasi-linear viscoelastic model has been introduced, 
which enables the experimental observations to be reproduced in a finite element 
calculation. For this purpose an internal variable has been used, which depends 
through its evolution equation on the deformation history. 

This modification violates the assumption of time invariance of the material 
response and is therefore inconsistent with the use of hereditary integrals (such 
as in the present quasi-linear viscoelastic model). Thus, the proposed modifi- 
cation cannot be considered as a valid extension of the quasi-linear viscoelastic 
approach. Rather it helps demonstrating the necessity of history dependent in- 
ternal variables for the description of the mechanical response of soft tissues in 
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repeated cycles. Non-linear hyperelastic visco-plastic models, such as the one 
proposed by Rubin and Bodner ([11]), should be used for this purpose. 

In the near future cyclic experiments on soft biological tissues will be per- 
formed in-vivo (before extraction) and ex- vivo (after extraction) , with particular 
attention to the differences in the softening attitude and changes in the internal 
structure of the tissue. 
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Abstract. Mechanical properties of biological tissues are needed for accurate 
surgical simulation and diagnostic purposes. These properties change post- 
mortem due to alterations in both the environmental and physical conditions of 
the tissue. Despite these known changes, the majority of existing data have been 
acquired ex vivo due to ease of testing. This study seeks to quantify the effects 
of testing conditions on the measurements obtained when testing the same tis- 
sue in the same locations with two different instruments over time. We will dis- 
cuss measurements made with indentation probes on whole porcine livers in 
vivo, ex vivo with a perfusion system that maintains temperature, hydration, and 
physiologic pressure, ex vivo unperfused, and untreated excised lobes. The data 
show >50% differences in steady state stiffness between tissues in vivo and un- 
perfused, but only 17% differences between in vivo and perfused tests. Varia- 
tions also exist in the time-domain and frequency domain responses between all 
test conditions. 



1 Introduction/Motivation: 

The mechanical properties of biological tissues are necessary for accurate surgical 
simulation and diagnostic purposes. Software-based simulation of surgical procedures 
relies on accurate representation of the mechanical response of tissues subject to sur- 
gical manipulations. If the tissue models or parameters are significantly different from 
reality, then negative training transfer may result from the use of the flawed simulator. 
Palpation, or manual evaluation of tissue mechanical response, has been used since 
the dawn of medicine, and numerous types of pathology cause changes in the viscoe- 
lastic character of tissue [1,2]. Should mechanical testing be used on this basis for di- 
agnostic purposes, accurate measurements of healthy and diseased tissues are needed. 

Soft tissues not only have complex material properties that are difficult to charac- 
terize, but also exist in an environment that affects their intrinsic behavior. Testing the 
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tissue in its natural state is ideal for ensuring accurate representations of the mechani- 
cal behavior we wish to characterize but difficult to achieve due to accessibility, ethi- 
cal, variability, noise, and uncontrolled boundary condition issues. Despite these is- 
sues, several groups have recently developed instrumentation to measure tissue 
properties in vivo [3—8] . Although these groups have successfully made force- 
displacement measurements on tissues in vivo, the interpretation of these results re- 
mains to be understood. 

The majority of existing data has been acquired under ex vivo conditions because 
this allows for precise control of boundary and loading conditions, provides access to 
appropriate testing sites and uses fewer animals [9-12]. However, testing soft tissues 
ex vivo drastically alters their properties and behavior [7, 13] and transplant research- 
ers indicate that tissues lose their functional viability and structural integrity within 
hours [14, 15]. A qualitative and quantitative understanding of the differences in 
measuring material properties from these different conditions is clearly needed. 

We believe that to understand the differences between testing conditions best, 
measurements that can capture both the elastic and viscous properties of soft tissues 
need to be made on the same organ, at the same location, across various environ- 
mental conditions. Several researchers have made soft tissue measurements in various 
environmental states [3, 5, 7] (in vivo, in situ, whole and partial organ ex vivo with 
and without controlling for hydration and temperature effects). Despite this, no exam- 
ples were found in which tissues were tested, harvested and retested in the same loca- 
tion with the same instruments to examine the changes in tissue properties post mor- 
tem. Brown et al [7] have measured the first squeeze force-displacement and stress 
relaxation response of porcine liver using graspers across three different environ- 
mental conditions: in vivo, in situ, and ex vivo. However, their data reveal structural 
not material properties with complex boundary conditions. Since they only measured 
the first squeeze response there is no measure of repeatability within location or 
across condition. 

This study seeks to quantify the effects of testing conditions on soft tissue material 
property measurements. We will discuss measurements made with two different in- 
dentation probes designed to capture the elastic and viscous material properties on 
whole porcine liver tissue in vivo, ex vivo with perfusion, ex vivo post-perfusion, and 
in vitro (i.e. warm ischemic partial organ tissue). These tests serve to verify the func- 
tion of our perfusion system, examine the differences between in vivo, perfused and 
unperfused tissue, and to compare intact versus excised organ conditions. 



2 Methods 

We seek to quantify the responses of tissues under four different conditions, namely 
the in vivo case, the in vitro excised lobe case, and two different ex vivo whole organ 
cases including perfused testing, and testing on tissues that have been flushed of 
blood with the perfusate, but tested thereafter without being supported by the perfu- 
sion system. 

The following sections describe the perfusion system in detail, as well as the inden- 
tation testing apparatuses used to measure the viscoelastic response of the tissues un- 
der each test condition. 
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2.1 Normothermic Extracorporeal Liver Perfusion System 



To accurately measure the mechanical properties of the liver, it is crucial that we 
maintain cellular integrity while keeping the organ in as natural a state as possible ex 
vivo. Thus we have built an apparatus similar in concept to normothermic extracopo- 
real perfusion systems using heparinzed Lactated Ringer’s solution as the perfusate 
(see Fig. 1). The system stores this solution in reservoirs suspended at specified 
heights to obtain the appropriate physiologic pressures into the hepatic artery (100- 
120mmHg) and portal vein (15-20mmHg). Both pressures and flow rates can be eas- 
ily adjusted by altering the height of the reservoirs and by partially closing tubing 
clamps respectively. The perfusate is then allowed to drain via the intrahepatic vena 
cava into a bath where it is heated to a physiologic temperature (39C for pigs) and 
circulated to the reservoirs via a pump. The solution also flows over the organ to 
maintain hydration without having to submerge the organ. To ensure consistency in 
our measurements, the organ rests on a sturdy plate covered with fine grit sandpaper 
to localize and stabilize the area of tissue under study, and the perfusion pressure is 
held constant rather than mimicking physiologic pulsatile pressure. 





Arterial reservoir 

Portal venous 
reservoir 

Data Aquitision 
system 



Porcine iver 



Fig. 1. (left) Our Normothermic Extracorporeal Liver Perfusion system schematic, (right) 
NELP system in use 



2.2 Test Instruments 

We used two different indentation devices to measure the viscoelastic response of the 
tissue. The TeMPeST (Tissue Material Property Sampling Tool) allows us to measure 
the small strain frequency response of tissues, while the VESPI (ViscoElastic Soft- 
tissue Property Indentation instrument) examines the large strain time-domain re- 
sponse (see Fig. 2). 
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Fig. 2. The TeMPeST and the VESPI devices testing 2 separate locations on the same pig liver 
maintained in the normothermic extracorporeal perfusion system. 



TeMPeST 1-D. TeMPeST 1-D is a 12mm diameter minimally invasive instrument, 
designed to measure the compliance of solid organ tissues within the linear regime. A 
5mm right circular punch vibrates the tissue while recording applied load and relative 
displacement. Mechanical bandwidth is approximately 80Hz when in contact with or- 
gan tissues, however the force and position sampling frequency is up to 2kHz, so 
measurements can be made to approximately 200Hz, depending on the material. 
Range of motion is 1mm and forces up to 300mN can be exerted. It has previously 
measured the properties of porcine liver and spleen in vivo , rodent (rat) liver and kid- 
ney ex vivo [6, 16], and has been used in initial investigations of bovine, ovine and 
human vocal tissue samples ex vivo. 



VESPI. The VESPI also performs normal indentation on tissues. It was designed for 
bench-top use, and subsequently modified to permit open surgical in vivo measure- 
ments as well. A 6mm diameter flat punch rests, with only 3g load due to counter 
weights, on the tissue surface until a standard laboratory mass is released (at zero ve- 
locity) onto a platform mounted co-axially with the indenter tip. The organ rests on 
the same platen to which the measurement arm and indenter tip are mounted. Loads 
from 20 to lOOg (nominal stresses of 6.2 to 31kPa) generate much larger deformations 
than those created by the TeMPeST. Because of the large range of motion of the in- 
strument, and the relative immobility of the organ resting on the platen (compared 
with TeMPeST testing), breathing does not need to be suspended during testing, ena- 
bling much longer data acquisition periods (typically 300 seconds). During this pe- 
riod, the angular position of the measurement arm is measured at a rate of 1kHz using 
a miniature contactless rotary position sensor (Midori America Corporation, CA) 
(resolution 1 3.7pm, signal to noise -100:1). Since the lever arm has a length of 
1 1.5cm, and the maximum depth of indentation was on the order of 10mm, small an- 
gle approximation was assumed and thus the voltage from the rotary sensor is con- 
verted directly to indentation depth using a linear gain. 
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2.3 Test Protocol 

In vivo tests were performed on deeply anesthetized animals on assisted ventilation 
with 100% oxygen. The TeMPeST instrument was used to acquire compliance data 
on either one or two locations on the liver during data acquisition periods of approxi- 
mately 20 seconds. During this time, ventilation was suspended to prevent pulmonary 
motions from saturating the position sensor measurements. Indentations using the 
VESPI device were made at the same location(s), but without the necessity for sus- 
pending ventilation. Organ thickness measurements were taken prior to every VESPI 
measurement with a 0-25 mm dial indicator. Initial position senor values were noted 
and loads were applied for 300 seconds. Once the load was removed the organ was al- 
lowed to recover to its preindented state (typically 200s, as was determined by com- 
paring the current position to the preloaded value). 

Following in vivo testing, heparin was injected systemically to prevent clotting, 
and the animal was sacrificed. The liver was harvested, and a lobe was removed and 
tested immediately with the TeMPeST (in vitro testing). The cut surface of the re- 
mainder of the organ was cauterized to prevent leakage and the organ was flushed 
with heparinized lactated Ringer’ s solution, packed in ice and transported to the labo- 
ratory. 

Upon arrival (60-80min post-sacrifice), the liver was connected to the arterial and 
hepatic venous perfusate reservoirs, and allowed to come to physiological temperature 
and pressure before testing was resumed. TeMPeST and VESPI tests were performed 
in the same locations as the in vivo tests to minimize variation in the measurements 
due to the unknown locations of large vessels or connective tissues within the organ. 
Testing on the excised (untreated) lobe was also performed over time with both in- 
struments. Times of sacrifice, and initiation and termination of perfusion were re- 
corded to permit examination of measurements over time. 

Following the completion of testing on the perfused organ (typically 2 hours), per- 
fusate flow was stopped, and the organ was tested again over time to observe any fur- 
ther changes in the response (typically 1 hour). 

The research was conducted in compliance with the Animal Welfare Act Regula- 
tions and other Federal statutes relating to animals and experiments involving animals 
and adheres to the principles set forth in the Guide for Care and Use of Laboratory 
Animals, National Research Council, 1996. 



3 Results 

Variations in the measured responses are observed between all of the conditions under 
consideration, including both changes in the measured tissue stiffness and the time 
dependent viscous character of the responses. The measurements performed with the 
TeMPeST and VESPI will be shown in the following sections. 

3.1 TeMPeST Results 

Fig. 3 shows the frequency response calculated from the ratio of the FFTs of the posi- 
tion and force signals, together with the ideal first order filter response of a Voigt tis- 
sue model. As the response appears to have a -20dB/decade slope after the break, and 
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the phase lag at high frequencies is approximately -90°, the Voigt/first order response 
is a reasonable first approximation to fit the results. It is observed that the in vivo 
measurements show the highest compliance (lowest stiffness), while the perfused tis- 
sues are stiffer, and the unperfused tissues (after prior perfusion) are the stiffest. In 
addition, the break of the response shifts to higher frequencies in each of these cases. 
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Porcine liver 5 location 1 compliance: perfusion effects 



Fig. 3. Frequency response of tissue measured with TeMPeST and Voigt model approximation 
of tissue response. First order filter characteristics include asymptotes to better show character- 
istic frequency (dashed lines). 



3.2 VESPI Results 

Fig. 4 shows representative results for three conditions made on the same liver taken 
at the same location. Strain is calculated as depth of indentation normalized with re- 
spect to thickness measured prior to that specific test. It can be seen that the tissue in 
vivo is softer than perfused tissue, which in turn is softer than the unperfused organ. 
Not shown is the data from the in vitro experiment, which was much softer than all 
other conditions and which never fully achieved a steady strain state. Most signifi- 
cantly, the perfused steady state strain is within 17% of the in vivo value, much closer 
than the unperfused strain, which differs by more than 50% (considering 3 rd indenta- 
tion for each case). 

We were also interested in examining the repeatability of the measurements within 
location. Fig. 4 shows that after allowing the tissue to recover fully after testing, the in 
vivo time responses are very similar to each other, as are the perfused tests. The un- 
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Fig. 4. VESPI response for three of the conditions under consideration on a 27kg pig liver at 
the same location. (A) In vivo lOOg load. The second and third indentations were taken 10 and 
50 minutes after the first. (B) Perfused lOOg load. The second and third indentations were taken 
20 and 48 minutes after the first. (C) Unperfused lOOg load. The first indentation was taken 1 
minute post perfusion while the second and third were 12 and 25 minutes after that. (D) Results 
of the first indentation performed at each condition. 



perfused test, because an external source of pressurized perfusate is no longer avail- 
able, does not recover, and significant changes are observed between the first and 
third measurements. 

Lastly, the curvature of the results provides insight to the creep time constants in 
the tissue. A quantitative analysis of the results will need to be performed, but qualita- 
tively it can be seen that the curves of the various conditions are indeed different sug- 
gesting that the testing conditions alter the viscous properties of the tissue. 



4 Conclusions and Future Work 

Humans have the ability to distinguish between seven different levels of unidimen- 
sional stimuli [17]. Despite this rather coarse just noticeable difference, it has also 
been shown that haptically humans can feel differences in force as small as 10% [18]. 
For this reason, employing measurements from unperfused tissue as a proxy for per- 
fused (or in vivo ) data may result in significant inaccuracies in surgical simulators for 
training. For example, if one becomes accustomed to manipulating virtual tissues that 
are stiffer than real ones, excessive, possibly damaging forces may be applied when 
the trainee reaches the operating room. We hypothesized that in creating an environ- 
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ment that closely approximates in vivo conditions, we can maintain the mechanical 
viability of the tissue such that the properties we measure ex vivo are comparable to 
the in vivo properties and that without such an environment, the material properties 
are altered far beyond 10%. 

Each instrument shows variations in the stiffness and damping properties of the 
liver depending on the test condition. However, the large deformation time responses, 
when using the perfusion apparatus, approach those of tissues tested in vivo. The pres- 
surized flow of perfusate permits the organ to fully recover after testing, just as the 
flow of blood through the living organ would. In addition the static pressure applied 
to the perfusate provides an internal “boundary condition” to the tissue which is pre- 
sent (at least on average) in the living organ, but is absent when testing tissues on the 
lab bench, even if immersed in physiological solution. Lastly, testing whole organs 
versus cut specimen provides a more accurate reference state containing residual 
stresses and strains. 

The TeMPeST measurements show the high frequency response of the tissue, 
showing a slight increase in the break frequency after the tissues have been perfused. 
One possible explanation is that the Lactated Ringer’s solution, mostly water, behaves 
as a Newtonian fluid, with a viscosity lower than that of the non-Newtonian blood 
that normally perfuses the tissues. Whether this is a noticeable error from the perspec- 
tive of simulation is a question that cannot be answered in this study. If diagnostic ap- 
plications are envisioned, it may be necessary to refine the perfusion system further to 
ensure a closer match in behavior over a wide range of time scales. 

More recent examination of the literature as well as input from a pathologist exam- 
ining some samples taken of the liver over time in each condition have shown that the 
perfusion pressure for the hepatic portal branch of the system is currently too high. In 
particular, the 20mmHg value is more than double the published value of 9mmHg 
[19]. This over pressure is most likely the reason why the shape of the VESPI creep 
curves between the in vivo and ex vivo perfused states were different. The increased 
pressure probably also accounted for some of the cellular dissociation seen histologi- 
cally. Measurements of in vivo portal venous and arterial pressure will be performed 
in future experiments and the static pressures of the system will be altered accord- 
ingly. Osmotic and oncotic pressure issues may also be playing a role. Discussions 
with transplant surgeons are underway to better improve our perfusate recipe for fu- 
ture tests. 

Tests scheduled, but not yet performed include testing whole organs over time in 
both the untreated (directly harvested) case and in the flushed-unperfused case. Also, 
a test will be conducted using the VESPI to determine the amount of mechanical 
damage (from analyzing histological samples) on the organ after the first load as a 
function of load. 

In conclusion, it was confirmed that untreated tissues behave much differently than 
tissues in vivo , and that the perfusion system provides a suitable environment for test- 
ing whole organ tissues. The maintenance of hydration, temperature, osmotic balance 
and internal pressure are necessary to permit extended testing of whole organs on the 
lab bench, where testing can be conducted more conveniently than in the operating 
room. In addition, the perfusion system will enable the use of organs harvested from 
other sources, without ever performing dedicated in vivo tests, reducing the cost of 
testing and the ethical and administrative issues of in vivo testing. It is recommended 
that researchers performing soft tissue testing on solid organs make use of similar sys- 
tems to provide mechanical function support to tissues tested in the laboratory setting. 
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Abstract. A Finite Element model of the face soft tissue is proposed 
to simulate the morphological outcomes of maxillofacial surgery. Three 
modelling options are implemented: a linear elastic model with small 
and large deformation hypothesis, and an hyperelastic Mooney-Rivlin 
model. An evaluation procedure based on a qualitative and quantitative 
comparison of the simulations with a post-operative CT scan is detailed. 

It is then applied to one clinical case to evaluate the differences between 
the three models, and with the actual patient morphology. First results 
shows in particular that for a “simple” clinical procedure where stress is 
less than 20%, a linear model seams sufficient for a correct modelling. 

1 Introduction 

Modeling the human soft tissue is of growing interest in medical and computer 
science fields, with a wide range of applications such as physiological analysis, 
surgery planning, or interactive simulation for training purpose [1]. In maxillo- 
facial surgery, the correction of face dismorplrosis is addressed by surgical repo- 
sitioning of bone segments (e.g. the mandible, maxilla or zygomatic bone). A 
model of the patient face to simulate the morphological modifications following 
bone repositioning could greatly improve the planning of the intervention, for 
both the surgeon and the patient. 

Different models were proposed in the literature. After testing discrete mass- 
springs structures [2], most of the authors used the Finite Element method to 
resolve the mechanical equations describing the soft tissue behavior. [3], [4] and 
[5] first developed linear elastic models. With more complex models, [6] discussed 
the advantages of non-linear hypotheses, and [7] began accounting for tissue 
growth in their simulation. 
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One of the most important issue in soft tissue modeling is to assess the quality 
of the simulations. From a modeling point of view, it enables to evaluate and 
compare different methods, for example linear versus non-linear models. This is 
above all essential for the surgeon since the use of a soft tissue model in actual 
surgical practice cannot be considered without an extensive clinical validation. 
While many models were proposed in the literature, few works propose satisfying 
validation procedures. 

In this paper, we first propose different modeling hypotheses of the face 
soft tissue. An evaluation procedure is then detailed, based on a qualitative 
and quantitative comparison of the simulations with a post-operative CT scan. 
Results are presented for a clinical case of retro-mandibular correction. The bone 
repositioning actually realized during the intervention is measured with accuracy, 
then simulated using the biomechanical model. Simulations are thus compared 
to assess the influence of modeling options, and their relevancy with respect to 
the real post-operative aspect of the patient. 

2 Modeling the Face Soft Tissue 

A project for computer-aided maxillofacial surgery has been developed for several 
years in the TIMC laboratory (Grenoble, France), in collaboration with the Pur- 
pan Hospital of Toulouse, France. In that context, a Finite Element model of the 
face soft tissue has been developed to simulate the morphological modifications 
resulting from bones repositioning. In [8] , we mainly presented our methodology 
to generate patient-specific Finite Element models. A generic mesh was built, 
organized in two layers of hexahedrons and wedges elements. The principle was 
then to conform this generic model to the morphology of each patient. Using 
elastic registration, nodes of the mesh were non-rigidly displaced to fit the skin 
and skull surfaces of the patient reconstructed from a pre-operative CT scan. 

Once a mesh of the patient is available, biomechanical hypothesis must be 
chosen to model the mechanical behavior of the face tissue. Three different meth- 
ods are compared in this paper: a linear elastic model, under small then large 
deformation hypothesis, and an hyperelastic model. 

2.1 Linear Elastic Model 

A first hypothesis is to model the tissue as a homogeneous, linear elastic material. 
This assumption, which considers the stress /strain relationship of the system 
as always linear during the simulation, is called mechanical linearity. Although 
biological tissues are much more complex, this behavior was found coherent for 
a relative strain under 10 to 15% [9]. The material properties can be described 
using the Hooke’s law with two rheological parameters, the Young modulus and 
the Poisson ratio. 

A second option of modeling depends on the deformations range. In small 
deformations hypothesis (also named geometrical linearity), the Green-Lagrange 
formula linking the stress and strain tensors is linearized by neglecting the sec- 
ond order term [10]. As a consequence, the formulation can be written as a linear 
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matrix inversion problem, which is straightforward and fast to resolve. Under 
the large deformations hypothesis, the second order term is not neglected, which 
leads to a more accurate approximation but dramatically increases the compu- 
tation complexity. 

A linear material with small deformations is the most widely used hypoth- 
esis in the literature of facial tissue deformation modeling. Despite the fact it 
is probably limited due to the complexity of the tissue properties and the sim- 
ulated surgical procedures, this model is certainly the first to be tested and 
compared with actual data. We had therefore implemented such a model, with 
rheological parameters of 15 kPa for the Young modulus, and a 0.49 Poisson 
ratio (quasi-incompressibility). Simulations were carried out under both small 
and large deformations hypotheses. 

2.2 Hyperelastic Model 

As stated in different papers, linear models become inaccurate when the dis- 
placements or the deformations are large [11], [6], and when the rheology of the 
material is non-linear. Numerical errors appears due to the non-invariance in 
rotation. A major shortcoming especially lies on the material constitutive law. 
Experiments on biological tissue [9] have shown that the stress increase much 
faster than the strain as soon as the small deformation context is not applicable. 
This increase of the stiffness must be taken into account, which is not possible 
with a linear constitutive law such as the Hooke’s law. 

Therefore, a classical modeling framework, the hyperelasticity, can be used 
to directly account for all the non linearities (mechanical and geometrical) in 
the mathematical formulation. Whereas a material is said to be elastic when 
the stress 5 at a point X depends on the values of the deformation gradient F, 
the material is said to be hyperelastic when the stress can be derived from the 
deformation gradient and from a stored strain energy function W : 

o dW 
d E 

where E is the Lagrangian strain tensor. 

The strain energy W is a function of multidimensional interactions described 
by the nine components of F. It is very difficult to perform experiments to de- 
termine these interactions for any particular elastic material. Therefore, various 
assumptions have been made to derive simplified and realistic strain energy func- 
tions. One of this assumption is the Mooney- Rivlin materials modelling [12]. For 
exemple, the energy function W can be approximated by a 5 coefficients Mooney- 
Rivlin material, so that: 

W = aio(/i — 3) + a2o(U — 3)“ + aoi(/2 — 3) + 002(^2 — 3) 2 + an(/i — 3)(/2 — 3) 

where I\ and I 2 are the first and the second invariant of the deformation tensor 
E. Assuming a constitutive law for facial tissues that is close to the constitutive 
law proposed by [13] for the human tongue, a two parameters Mooney-Rivlin 
material was finally assumed for the simulations : aio = 2500 Pa and 020 = 625 
Pa. The three other parameters Oqi, oq 2 and an are set to zero. 
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3 Validation Procedure 

Few authors have proposed extended validation procedures for soft tissue model- 
ing. In maxillofacial surgery, most of them compare their simulations with facial 
and profile pictures of the patient. While a qualitative comparison is always re- 
quired, this method is quite inaccurate and does not afford a real tri-dimensional 
evaluation. The main other approach rely on the acquisition of the post-surgical 
patient morphology with an optical laser scanner. This enables a 3D quantitative 
comparison [4]. However, it is very sensitive to the accuracy of the acquisition 
and to the registration procedure that expresses it in the pre-operative patient 
referential. Moreover, there is always an important error between the simulated 
intervention and the bone repositioning actually realized during the surgery. 
The most advanced quantitative evaluation was recently proposed by [7], who 
measure the distances between their simulations and a post-operative CT scan. 

The evaluation protocol we propose also requires the acquisition of a pre and 
a post-operative CT scan. While a post-operative exam is invasive in terms of 
radiations, it is clearly the best available data to assess the quality of numerical 
simulations. With the improvement of modern scanners, its use can therefore be 
acceptable in a research context. 

Our evaluation procedure consists in four steps: 

1. measuring the bone repositioning actually realized during the surgery, by 
direct comparison of the pre- and post-operative data; 

2. simulating the bone osteotomies and applying the measured displacements 
to the bone segments; 

3. simulating the resulting soft tissue deformation using the biomechanical 
model; 

4. evaluating the differences between the simulation and the post-operative 
data, both qualitatively and quantitatively. 

The first two steps are realized using mathematical tools initially developed 
for a 3D cephalometry project [14] . Anatomical landmarks in areas that are not 
modified during the surgery are defined in both the pre- and post-operative CT 
slices, to register the two datasets in a same referential. Then, landmarks lo- 
cated on each bone segment (e.g. the mandible and maxillae) enable to measure 
the displacements actually applied during the surgery (figure 1). Although the 
anatomical landmarks are manually positioned on the CT slices, it has been 
shown that their repeatability is in mean .25 mm, which yields to very accept- 
able results in the measurements of the bone displacements. Moreover, a rigid 
registration can be added to further improve the accuracy the measurements. 

The measured displacements define the boundary conditions for the Finite 
Element model. Inner nodes in contact with the non-modified skeleton surface 
are fixed, while the measured displacements are applied to the nodes on the os- 
teotomized bone segments. Nodes around the osteotomy line are not constrained, 
to account for the bone-tissue separation of the surgical access. Rest of the nodes, 
in the outer part of the mesh and the mouth and cheeks area, are let free to move. 
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Fig. 1. Clinical case of mandibular prognatism. Left, the patient skull and skin surface, 
before and after the surgery. By comparison of the two skeleton surfaces using our 3D 
cephalometry, the mandibular correction actually realized during the intervention was 
accurately measured and reproduced (right). 



Once the outcome of the surgery has been simulated with the biomechanical 
model, it can be compared with the post-operative skin surface of the patient, 
reconstructed from the CT scan. The quantitative comparison between the two 
datasets is achieved using the MESH software [15], which has been improved to 
calculate signed Euclidian distances. 

4 Results 

Results are presented on a clinical case of retro-mandibular correction. A pre- 
operative CT scan was first acquired, which enabled us to generate a 3D mesh 
conformed to the patient morphology. After the patient was operated in a conven- 
tional way, a post-operative CT scan was acquired. By comparing both datasets, 
the actual displacement applied to the mandible during the intervention was 
measured (figure 1). It consisted in a backward translation of 0.9 mm (in the 
mandible axis), and a slight rotation in the axial plane. The measured proce- 
dure was then reproduced on the skeleton model, and boundary conditions for 
the Finite Element model were set. 

By comparing the vertebras positions in both CT scans, it can be seen that 
the inclination of the head was not similar during both exams, with a difference 
of more than 10 degrees. Unfortunately, this imply the simulations would not 
be comparable with the post-operative aspect in the neck area, since the head 
position directly influence the cervico-mentale angle. Therefore, other boundary 
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conditions were added to the extreme lower nodes of the mesh, in the neck area, 
to reproduce this modification of the head position. Although quite qualitative, 
this should enable us to better compare the simulations with the actual patient 
morphology. 

Simulations were computed using the Ansys™ Finite Element software (An- 
sys Inc.). For the linear elastic model, the computing time is less than 3 seconds 
with the small deformation hypothesis, and almost 3 minutes in large defor- 
mations. The hyperelastic calculus required up to 8 minutes. All simulations 
were static, and ran on a 2.4 GHz PC. The model counts around 5000 elements 
and 7650 nodes. For the two non-linear models, numerical convergence can be 
difficult to obtain if the boundary conditions are not well defined. This was par- 
ticularly the case in our first trials, before the mesh was extended in the posterior 
direction, which enabled us to better constraint the posterior nodes. Generally 
speaking, it must be recall that convergence is very sensitive to the boundary 
conditions, the quality of the elements shape and the time-steps used during the 
numerical resolution. 

Figure 2 shows the Von-Mises repartition of the strain, calculated for linear 
model in small deformation. Figure 3 presents the results obtained with the linear 
elastic model in large deformation, along with the post-operative skin surface of 
the patient. Finally, figure 4 shows the distances measured between the models 
predictions and the patient data with the MESH software. 



5 Discussion 

Before analyzing the results, a methodological point should be discussed. The 
use of numerical data (CT scan) enable us to obtain a quantitative evaluation 
of the simulation errors. Such results are quite rare (only [7] got similar ones 
with actual post-operative data) and seems extremely important and necessary 
to really assess the influences of the different modeling options. Nevertheless, the 
numerical values should be carefully analyzed. They represent the minimal Eu- 
clidian distance between points of the model and the patient skin surface, which 
does not mean the distances to their real corresponding points (e.g. the distance 
between a point P and the surface S can be smaller than the distance between 
P and its actual corresponding point P' in S ) . These numerical values are thus 
always a minimization of the true errors. Therefore, the quantitative analysis 
must always be completed by a qualitative evaluation of the results, carried out 
by a clinician. This remains the best way to know how well the model is per- 
ceived by the surgeon, an gives an emphasis to the most relevant morphological 
areas in the face: cheeks bones, lips area, chin and mandible angles. 

First, it should be noted that the simulations obtained with all models are 
quite similar. This is an interesting result since the strain repartition (figure 2) 
shows that the relative strain are above 20% in a large area around the mandible, 
with peak of 30% to 100% in the osteotomy area (the maximum being in the 
bone surface region). A first conclusion is thus that a linear model in a small 
deformation framework appears quite acceptable even for relative deformation 
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Fig. 2. Repartition of the strain for the linear model in small deformation. 




Fig. 3. Comparison of the post-operative data (left) with the linear elastic model in 
large deformation (center). Both models are superposed in the right view. 




Fig. 4. Measurement of the error between the simulations and the post-operative data, 
for the linear model in small (left) and large deformation (center), and the hyperelastic 
model (right). 
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up to 20%. The large deformation hypothesis does really not decrease the errors. 
Surprisingly, the first results obtained with the hyperelastic model shows more 
important errors. This could be explained by the fact this modeling is much 
more sensitive to several critera like the mesh (which must be refined in the 
high-stress areas), the boundary conditions and the rheological parameters. Such 
complicated models require more testing before being used, and may not be the 
most adapted for problems with relatively small deformations. 

Clinically, the simulations are of good quality and quite coherent with the 
actual outcome of the surgery. The accuracy is the best in the chin area, which 
seems logical since that region is one of the most constrained. A slight swelling 
is observed in the cheeks area of the model, which is a known clinical behavior 
in retro-mandibular procedures that was correctly reproduced with the model. 

Although they are not the largest numerically, between .5 and 2 mm, errors 
around the lips are important. It can be observed that the shape of the infe- 
rior lips is unchanged from the pre-operative state, just translated, and thus 
incorrect. This is explained by the fact contacts between both lips and between 
the lips and the teeth are not taken into account so far. Indeed, penetration 
occurred with the teeth, which would certainly have modified the shape of the 
lip if the contacts were handled. This essential modeling aspect, not discussed 
in the literature of facial simulation, will really have to be integrated. 

The most important errors, up to 6 mm, occur around the angles of the 
mandible. Numerical values should be analyzed carefully since the difference of 
head inclination certainly influence the results. Nevertheless, errors are expected 
important in that areas where the stress is maximum. They correspond to the 
osteotomy region, and thus the frontier between constrained and free nodes. 
The swelling observed in all models is more important that in the actual data. 
Before complicating the model, for example with growth modeling [7], we prefer 
to continue the tests and evaluations, since the behavior in these areas appears 
quite sensitive to the boundary conditions, especially for the two non-linear 
models. 

6 Conclusion 

A qualitative and quantitative evaluation procedure was proposed in this paper 
and used to compare different modeling of the face soft tissue. First results are 
quite encouraging. It has particularly been shown that for maxillofacial surgery 
simulation, a linear elastic model can be sufficient for simple procedures like 
retro-mandibular correction. 

Future works are to extend the evaluation of different modeling options and 
to assess the influence of elements (refinement, linear or quadratics elements, 
etc.), rheological properties and numerical methods. Lips and lips-teetlr contact 
must also be taken into account. Two more complex clinical cases are planned 
for the evaluation (with a post-operative CT scan): a bimaxillary correction with 
genioplasty, and a distraction of the orbito-zygomatic structure. The non-linear 
models are expected to be necessary to simulate these difficult, large deforma- 
tions procedures. 
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Abstract. The biomechanical properties of soft tissue derived from experimen- 
tal measurements are critical for developing a reality-based model for mini- 
mally invasive surgical training and simulation. In our research, we focus on 
developing a biomechanical model of the liver under large tissue deformation. 
This paper presents the experimental apparatus, experimental data, and 
formulations to model the experimental data through finite element simulation 
and also compare it with the hyperelastic models in the literature. We used 
tissue indentation equipment to characterize the biomechanical properties of the 
liver and compared the local effective elastic modulus (LEM) derived from 
experimental data with that from plane stress and plane strain analysis in 
ABAQUS. Our results show that the experimentally derived LEM matches 
closely with that derived from ABAQUS in plane stress and plane strain 
analysis and the Ogden hyperelastic model for soft tissue. 



1 Introduction 

In surgery, probing soft tissue is one of the most common tasks to ascertain the tissue 
characteristic as being hard or soft. Hence, reality-based modeling of soft tissues is 
critical for providing accurate haptic feedback to the surgeon in surgical training and 
simulation. By reality-based modeling, we are interested in modeling tissues as accu- 
rately as possible by determining the mechanical properties experimentally. In the 
literature, most modeling efforts assume the mechanical properties of the soft tissue 
and develop methods to efficiently solve the tissue simulation problem for robot- 
assisted surgery/training. However, the experimentally measured material properties 
are critical for accurate tissue modeling. The goal of this paper is to derive the local 
effective modulus (LEM) of the liver tissue as it is compressed over a large range and 
compare the experimental results with the finite element simulation in ABAQUS and 
some hyperelastic models in the literature. In our study, pig liver was used as the 
sample tissue for deriving the material properties. The technique developed in this 
paper can be easily extended to characterize the material properties of other soft tis- 
sues as well. 
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“Global” elastic deformations of real and phantom tissues have been studied exten- 
sively in previous work, through simple poking interactions [1], However, these 
methods are simplistic since they do not consider the complex boundary conditions 
that are normally present, both internal to the organ and on the exterior surface. Howe 
and colleagues [2] have developed a “truth cube” for validation of models, however 
they have not yet extended this model to tool-tissue interactions for common surgical 
tasks such as probing tissues. There has also been research on estimating the me- 
chanical properties of the tissue through high-frequency shear deformations of the 
tissue sample, and elastography techniques. A variety of other techniques also exist in 
the literature for estimating the viscoelastic characterization of tissues, for example, 
[3, 4], Ottensmeyer [5] and others have performed tissue experiments to characterize 
force vs. displacement for pig liver tissue, however the tissue displacement was small 
and they have not characterized the mechanical properties of the liver. The theory of 
elastic material subjected to large deformations has been used to model the physical 
behavior of the biological tissues [6]. 

The quantitative knowledge of the biomechanical property of tissue is essential for 
soft tissue modeling. Fung [7] showed that the elasticity property of rabbits’ mesen- 
tery could be simply expressed as an exponential function. To develop a liver model 
for probing, it is necessary to characterize the material properties of the liver over a 
large range of deformation consistent with the range of deformation in surgery. 
Hence, it is necessary to derive the biomechanical properties of the liver, which is 
valid for both small and large strain regions [8]. The goal of this paper is to develop a 
computational model in ABAQUS and verify how well the experimental data fits 
with: a) ABAQUS computed stress vs. strain relationship for a specimen and b) some 
hyperelastic models in the literature. 

The rest of this paper is organized as follows. In section 2, we describe the materi- 
als and methods used to derive the LEM for the liver over a large deformation range. 
This is followed by the mathematical formulation of the large probe model and the 
various hyperelastic models that we will validate with the experimental data and 
ABAQUS model. We present the computational method in ABAQUS to derive the 
LEM over the range of deformation using both plane stress and plane strain analysis. 
In section 3, we present our experimental results and compare the various hyperelastic 
models with the experimental and ABAQUS computed material properties. Linally in 
section 4, we make some concluding remarks and discuss our future work in this area. 



2 Materials and Methods 

2.1 Experimental Setup 

We have designed and developed a tissue compression apparatus, which can measure 
the compressive force and displacement. Ligure 1 shows the configuration of our 
experimental system. The system consists of a motion control part, a force measuring 
part, and a post-data processing part. The motion control part is a lead screw assem- 
bled with a geared DC motor and encoder (Maxonmotor, Inc.), which is supported by 
two horizontal supports. The antibacklash nut in the lead screw prevents any backlash 
in the mechanism. A precision JR3 6 axis force/torque sensor (model 85M35A-I40) 
was attached to the probe and it travels along the lead screw as shown in Ligure 1. 
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When the motor rotates, the probe moves 
as per the motion control command. The 
position of the probe was controlled by 
the dSPACE DS1103 controller board 
(dSPACE, Inc.), which also records the 
force and displacement data. The sam- 
pling frequency for control and data 
acquisition was 1000 Hz. The algorithm 
implemented in the dSPACE card was a 
proportional +derivative (PD) control 
scheme. The communication between the 
dSPACE card and the JR3 data acquisi- 
tion card was accomplished by CLIB 
library (dSPACE, Inc.). The size of the 
probe was 50mm x 50 mm and the sur- 
face was polished and covered with pe- 
troleum jelly to minimize the contact 
friction force. The size of the probe en- 
sured that the liver sample fully con- 
tacted with the probe surface and contact 
over the entire surface was maintained as 
the tissue was compressed. We also assumed that over the entire range of compres- 
sion, the volume of the sample was conserved (though in reality this is not generally 
true). The liver samples were taken from freshly slaughtered pigs and transported to 
the lab within 2 hours post mortem. A cubic sample with the size of 10 mm x 10 mm 
x 10mm was cut from the liver. In the compression experiment, the liver sample was 
compressed to 30% of its nominal thickness, and the probe speed was 6.096mm/min. 
We selected the low velocity of the probe for quasi-static analysis of tissue deforma- 
tion. 




Fig. 1. Experimental setup for measuring 
force and displacement during tissue 
compression tests. 



2.2 Large Probe Model 

To determine the tissues characteristic as being hard or soft in surgery, the surgeon 
probes the tissue with a probe, which is significantly smaller than the organ surface. 
From an analysis viewpoint, small probe analysis is a nonlinear boundary-value prob- 
lem and the corresponding numerical problem is difficult to solve. One possible ap- 
proach for developing a numerical model for small probe analysis is to develop a 
working model from large probe analysis, which derives the stress and strain relation- 
ship from the force and displacement data. In the large probe analysis, the size of the 
probe is much larger than the sample, which ensures the uniform deformation of the 
sample and consequently simplified boundary condition. The material properties 
derived from the large probe analysis will eventually be used to develop a numerical 
model for small probe analysis. For our initial approach, we have assumed the tissue 
to be incompressible, homogeneous, and isotropic elastic material. There are several 
strain energy potentials that can be used to model soft-tissue deformation, such as the 
Fung [7] and Ogden model [6] etc. 
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2.2.1 Experimental Analysis 

Assuming the probe force as F, the 
stretch ratio as X, and the initial contact 
area of the cube as A 0 , the Cauchy stress 

F , 

O is given by: <7 = — A and the strain 

A 

l 




vide the experimentally observed force 
versus displacement plot into small sub- 
regions (see Figure 2), then we can as- 
sume that each subregion is linear. For 
each such sub-region, the stress is Oj and 




Fig. 2. Raw experimental data from liver 
probing. 



the strain is e :. Based on the assumption 
of linearization, we hypothesize that there is an experimental LEM, which relates the 
incremental stress to incremental strain. This expression is given: 



A o= A s*E j => E i = 



(7, - <T 



£ - £ 



for i = 1,2 
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2.2.2 Hyperelastic Models 

For large probe analysis 
(see Figure 3), the bound- 
ary conditions are rela- 
tively simpler. The dis- 
placement of the tissue’s 
top surface is assumed to son Tissue 

be the same as the dis- 
placement of the probe. 

There is no lateral stress on 
the top and bottom sur- 
faces. The sides of the 

tissue are also assumed to be stress free. The Piola-Kirchhoff stress tensor is related to 
the stress tensor, S, at a given point by the expression: P = SF‘ T . Flence the Piola- 
Kirchhoff stress tensor is given by: 



Compressed State 
Large Probe 

ism 



77- 



Fig. 3. Schematic of tissue compression for large probe model. 
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where p is the Lagrangian multiplier corresponding to the incompressibility con- 
straint, Det(F) = 1 , B = FF 1 is the Cauchy-Green left dilation tensor, and 0) ] = — — 

dl„ 



and (O. = 



9U 

dlF 



where / = trace(B) and II B = 
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trace(B) -trace(B) 
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principal invariants of B. For an incompressible and isotropic material, the strain 
energy function, U is a function of the principal invariants of B, hence, U = \J(I B ,II B ). 
The Piola-Kirchhoff stress tensor, P, hence satisfies the equation: 

Div (P) = 0 (3) 

From equation (2) and (3), we see that, the constitutive equation for the stress can 
be determined if the strain energy function is known. Monney-Rivlin model and 
Ogden model are two famous strain energy functions for modeling the rubber-like 
materials. In the following part, we will derive the constitutive equation for the stress 
respectively from these two models. 

Mooney-Rivlin Model: The strain energy function for the Mooney-Rivlin Model 
(valid for rubber-like materials) is given by: 

U=C l0 (I B -l) + C 01 (II B -l) (4) 



where U is the strain energy per unit of reference volume and C 10 , and C 01 are tem- 
perature-dependent non-negative material parameters. Thus the stress vs. strain rela- 
tionship for the Mooney-Rivlin model is given by: 

0 = 2(1 -e 3e )(C 10 e- £ +C 0l ) (5) 



Ogden Model: The strain energy function for the Ogden model is given by: 



U 



= Z^(a;' + ^ , + ^ , -3) 

;=1 a ; 



(6) 



where U is the strain energy per unit of reference volume and |i ( , and a, are the mate- 
rial parameters. Thus the stress vs. strain relationship for the Ogden model is given 
by: 



<7 = 



2^-(e~ £0r,+e 
i=1 «,■ 



1 

-eoCi+e 



) 



(7) 



2.2.3 Finite Element Modeling (ABAQUS) 

It is desirable to construct a predictive computational model that can simulate the 
tissue probing process and predict the mechanical response (probing force versus 
displacement characteristics) of liver probing. We propose to use twodimensional 
finite element (FE) models for this purpose, a plane-stress FE model and a plane- 
strain FE model. 

A 2D finite element model was built in ABAQUS (Release 6.3-1) to simulate the 
compression process and iteratively obtain the computational LEM. Figure 2 shows 
the experimental force vs. displacement data for a cube sample undergoing compres- 
sion. In the analysis below, we subdivided each force vs. displacement curve into sub- 
region and each region was assumed to be linear. Figure 4 shows the corresponding 
model in ABAQUS during tissue compression. The cube was defined as the deform- 
able material and composed of 90 elements. The probe was defined as a rigid body 
and the friction on the contact surface was assumed to be zero. The central line of the 
model was constrained vertically so that the cube would not slip due to the remaining 
stress. 

Two element types, plane stress and plane strain, were tested separately to find the 
local effective modulus. In the 2D plane stress model, the element type was CPS4, a 
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4-node bilinear element and in the 2D 
plane strain model, the element type 
was CPE4, again a 4-node bilinear 
element. The displacement of the 
probe was used as the input and the 
reaction force of the probe was com- 
pared with the corresponding experi- 
mental force for a given tissue dis- 
placement. The total displacement (up 
to 30% deformation) was divided into 
60 sub-regions (each sub-region had 
0.5% strain increment). Similarly, we 
also did the analysis where the 30% 
deformation was divided into 30 sub-regions such that each sub-region corresponded 
to 1% strain increment. 

For the large deformation, the deformed geometry has an obvious effect on the 
strain, i.e., as the tissue in progressively deformed the history of deformation has an 
affect on the modulus at a given displacement of the tissue. Thus in the ABAQUS 
formulation, for any given sub-region j (j = 1,2, ...) of the force vs. displacement 
curve, the current geometry of the deformed mesh has to be imported into the model 
from the previous step to compute the values for the new linear region. This was ac- 
complished by using the IMPORT option in ABAQUS, which allows us to import the 
geometry and the mesh of the model. The initial stress for each increment was set to 
zero because the initial stress does not affect the LEM for the increment (due to plane 
stress and plane strain analysis). We assumed Poisson’s ratio of 0.3 and an initial 
LEM of arbitrary magnitude E t ■ to start the simulation. Then the experimentally 




Fig. 4. Compression test using large probe. 



measured displacement AU exp was applied to the node that models the large probe. 
The linear FE analysis is performed and the FEcomputed reaction force AF fem of that 
node was compared with the experimentally measured AF exp . The E value is updated 
by equation (8) until AF fem of the new iteration converged to the experimentally 
measured AF exp as per the convergence criterion given by Equation (9). The con- 
verged value for the local effective modulus was denoted as E effective for that region. 

k F expt 

= for i,j= 1,2,... (8) 

A F 



II AF fem - AF fem II 

" i ^ ^ 0-02 for .7=1,2,... (9) 

DF exp 

J 

The first subscript, i, of E t -, is the iteration number in each sub-region to compute 
LEM, and the second subscript, j, is the sub-region number. Once LEM is computed 
for a sub-region, j, this is the initial input value of LEM for the next sub-region. The 
deformed geometry and the meshes for the sub-region j are imported into the model 
for the subregion (/+1) and the process is repeated until the final sub-region is ana- 
lyzed. 

This iteration procedure is schematically shown in Figure 5. The local effective 
modulus E effective so determined is a measure of the deformation of the tissue consis- 
tent with the experimental data and the FE discretization. One application is to subse- 
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quently form an FEM mesh of the 
liver, assign the elements in the 
FEM model with their respective 
^effective, an£ J use t h at pgM model 
to virtually simulate liver probing 
and vary the probing parameters 
such as the probing speed. Such a 
FEM model embedded with self- 
consistent LEM would be able to 
predict the monotonic deformation 
segments of the probing force 
versus displacement characteristics 
consistent with experimentally- 
measured values, should actual 
experiment of that particular prob- 
ing be performed. 



3 Results 

To understand the effect of the 
strain increment in the computa- 
tion of the local effective modulus, 
we did the FEM simulation in 
ABAQUS using two different 
strain increments, namely, 0.5% 




Fig. 5. Flowchart for LEM computation. 



and 1%. With 0.5% strain increment sub-division, we have twice as much data for 
LEM compared to 1% strain increment. For clarity of the figures below, we have 
shown the computed values for intermittent strains in the specimen. 



3.1 Comparison of LEM Computed with 0.5% and 1% Strain Increment 

Figure 6(a) shows the representative plot of LEM computed with 0.5% and 1% strain 
increment at a given displacement of the probe. As seen from the plot, the LEM val- 




Fig. 6. a) Representative plot of the LEM computed with 0.5% and 1% strain increment. Local 
effective modulus computed from FEM with plane stress and plane strain analysis for: b) 0.5% 
increment in strain and c) 1% increment in strain in FEM simulation. 















Characterization of Soft-Tissue Material Properties: Large Deformation Analysis 35 



ues for 0.5% and 1% strain increments are very close to each other for a given dis- 
placement. Furthermore, Figures 6(b) and 6(c) shows the computed LEM value in 
both plane stress and plane strain mode for 0.5% and 1% strain increment. Since the 
tissue was compressed up to 30% of its nominal height, there were twice as many data 
points in 0.5% increment of strain compared to 1% increment. However for the sake 
of clarity, fewer data points were plotted in figures 6(a-c). We have also computed the 
errors between the experimentally computed LEM and that obtained from ABAQUS 
in both plane stress and plane strain modes and observed this error to be less than 
10%. This leads us to conclude that the two dimensional ABAQUS model can pro- 
vide a reasonable estimate of the local effective modulus in the tissue over the range 
of compression. Based on the observed closeness in the local effective modulus for 
0.5% and 1% strain increment from ABAQUS simulation, 1% strain increment for 
analyzing other tissue samples should be adequate. 



3.2 Variability of LEM across Specimens 

We computed the LEM 
values at different displace- 
ments of the soft tissue for 
four different specimens. 

Figure 7 shows the LEM 
values obtained from ex- 
perimental analysis and also 
from ABAQUS. In the fig- 
ure plane-stress values have 
been plotted, however the 
planestrain values were very 
similar to the plane-stress 
values. As seen from the 
figure, the values for LEM 
were within close range of 
each other across different 
tissue samples. Based on the 
data obtained from the 
ABAQUS analysis, we 
plotted the results for all the 
four samples (all done at 1% 
strain increment) in matlab. We did a polynomial fit in matlab to determine the rela- 
tionship between the LEM and the displacement of the tissue. The mathematical ex- 
pression for the relationship between LEM, E, and the displacement of the tissue, x, 
for samples 1 through 4 is given by: 

Sample 1 :E= 0.0054a 4 - 0.014a 3 + 0.025a: 2 - 0.033a + 0.0054 
Sample 2: E = 0.0016a 4 - 0.00082a 3 + 0.012a 2 - 0.016a + 0.019 
Sample 3: E = 0.0014c 4 - 0.0009a 3 + 0.0098a 2 - 0.012a + 0.014 
Sample 4: E = 0.00 16.x 4 - 0.051a 3 + 0.077a 2 - 0.027a + 0.0086 



a) 




c) 



E Modulus of sample 3 




b) 




d) 



E Modulus of sample 4 




Displacement (mm) 

Fig. 7. Experimental and ABAQUS LEM for four samples. 
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Fig. 8. Stress vs. strain curve fit for 
the Ogden and Mooney-Rivlin model 
with the ABAQUS data. Experimental 
stress vs. strain curve is also shown. 



Table 1 . Mooney-Rivlin and Ogden model parame- 
ters for different tissue samples. 





Parameter 


Sample 

1 


Sample 

2 


Sample 

3 


Sample 

4 


Mooney- 

Rivlin 

Model 


c 

^10 


0.039 


0.067 


0.052 


0.063 


r 

''"01 


-0.041 


-0.067 


-0.052 


-0.066 


Ogden 

Model 


fi 


-0.228 


-0.221 


-0.179 


-0.380 




0.162 


0.230 


0.177 


0.272 




0.067 


-0.002 


0.007 


0.109 


«i 


11.085 


7.625 


6.663 


13.514 


a 2 


11.265 


7.913 


7.062 


13.616 


a 3 


10.996 


-25.000 


5.846 


13.465 



3.3 Comparison of Hyperelastic Models with Experimental Analysis 
and ABAQUS Computations 

Finally, we were interested in comparing the various hyperelastic models with the 
experimentally observed and ABAQUS computed stress vs. strain curve for large 
tissue deformation. Figure 8 shows the plot for a particular sample. We did this analy- 
sis for four samples. Table 1 shows the various parameter values for the Mooney- 
Rivlin and Ogden model. As seen from the table, the second parameter for the 
Mooney- Rivlin model is negative for all four samples. The negative parameter value 
for C01 in the Mooney- Rivlin model is not physically valid, as the parameters of the 
Mooney-Rivlin model have to be non-negative. In addition, the Mooney-Rivlin model 
is only valid for rubber-like materials. On the contrary, the experimentally observed 
and ABAQUS computed stress vs. strain curve fits well with the Ogden model. How- 
ever, it should be noted that agreement of the Ogden model is valid only for qua- 
sistatic analysis presented in this paper. Validity of this model for normal probing 
speeds is the area of future work. 



4 Conclusions and Future Work 

We have designed and developed a tissue compression apparatus for characterizing 
the mechanical response of liver tissue using large probe analysis. In our initial work, 
the liver tissue is assumed as incompressible, isotropic, and homogeneous elastic 
material. Based on the results from the compression test, we have developed a large 
probe model for the liver tissue. The tissue was compressed up to 30% of its normal 
height. We did the analysis with 0.5% and 1% strain increment in ABAQUS and 
concluded that the 1% strain increment in the computations produced satisfactory 
results. We built a 2-D finite element model to simulate the tissue probing process and 
computed the mechanical response during liver probing. We conducted a simulation 
with both plane-stress and plane-strain analysis. The ABAQUS computed LEM of the 
liver tissue was iteratively obtained. The results showed that the experimental and 
computational LEM were very close to each other for a given displacement of the 
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tissue. We also computed the LEM values at different displacements of the soft tissue 
for four different specimens. The values of the LEM were within close range of each 
other across different tissue samples. We also performed a polynomial fit in Matlab to 
determine the relationship between the LEM and the displacement of the tissue. Fi- 
nally, we compared the various hyperelastic models with the experimentally observed 
and ABAQUS computed stress vs. strain curve for large tissue deformation. Based on 
the results of all four samples, we concluded that the experimentally observed and 
ABAQUS computed stress vs. strain curve fits well with the Ogden model. 

The work presented in this paper represents the initial step toward developing a re- 
ality-based model for tool-tissue interaction during probing. The results presented in 
this paper are valid only for quasi-static analysis. Our next step would be to compute 
the LEM over a range of probing speeds. The computational model developed with 
this research will be the basis for solving the complex small probe problem. 
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Abstract. Advancements in robotics have led to significant improvements in 
robot-assisted minimally invasive surgery. The use of these robotic systems has 
improved surgeon dexterity, reduced surgeon fatigue, and made remote surgical 
procedures possible. However, commercially available robotic surgical systems 
do not provide any haptic feedback to the surgeon. Just as palpation in open 
procedures helps the surgeon diagnose the tissue as normal or abnormal, it is 
necessary to provide force feedback to the surgeon in robot-assisted minimally 
invasive procedures. Therefore, a need exists to incorporate force feedback in 
laparoscopic tools for robot-assisted surgery. This paper describes our design of 
a laparoscopic grasper with tri-directional force measurement capability at the 
grasping jaws. The laparoscopic tool can measure grasping forces and lateral 
and longitudinal forces, such as those forces encountered in the probing of tis- 
sue. Initial testing of the prototype has shown its ability to accurately character- 
ize artificial tissue samples of varying stiffness. 



1 Introduction 

The introduction of robot-assisted surgery into the operating room has revolutionized 
the medical field. The use of these systems not only has the advantages of traditional 
minimally invasive surgery (MIS), such as reduced patient trauma and recovery time, 
lower morbidity, and lower health care costs, but they also eliminate surgeon tremor, 
reduce the effects of surgeon fatigue, and incorporate the ability to perform remote 
surgical procedures. However, current systems have shortcomings when compared to 
traditional MIS, such as high cost, inability to use qualitative information, and lack of 
haptic feedback to the surgeon [1], Thus, the loss of force feedback capability for the 
surgeon has motivated several researchers to explore possible methods of restoring 
this feature that is commonplace in open procedures and to a lesser extent in laparo- 
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scopic procedures. Several researchers in this field have proposed solutions to incor- 
porate force feedback into current laparoscopic tools [2-8], These solutions have 
involved adding force sensors, such as strain gages, to record the forces at the tool tip 
through indirect methods. Researchers have also added force feedback by creating 
new laparoscopic tools or systems with sensors included in the designs [9-14], These 
designs have placed sensors away from the tool tip and either indirectly measure the 
forces at the tip or only measure one resolved force on the tool. Researchers have also 
incorporated a direct sensing method for tissue characterization through pressure 
measurements normal to the surface of the jaws [15, 16]. However, these methods are 
expensive, non-sterilizable, and not modular, which make them difficult to incorpo- 
rate into laparoscopic tools. 

While there have been significant advances towards incorporating force feedback 
into laparoscopic tools, there still exist problems that must be considered. Backlash 
and “play” within mechanisms can lead to inaccurate force measurement and feed- 
back through incorrect positioning of the end effector of these tools. The lack of fric- 
tion modeling can also lead to inaccuracies; especially since commercial laparoscopic 
tools and many newly developed prototypes use indirect measurement techniques for 
measurement of end effector forces. While increasing the feedback gain by several 
times would overcome the friction in the mechanism, it would also lead to in-exact 
forces felt by the surgeon and, consequently, make the overall system less transpar- 
ent. One of the most common approaches in measuring forces during grasping is 
through placing the sensors on the tool shaft or at the handle. This is an indirect ap- 
proach of force measurement. Also, the forces measured in these approaches are a 
resultant force in the grasping jaws and not the individual components of forces in 
grasping tissues and/or organs. Indirect force measurement is performed through a 
calibration between the sensor or driving motor current and the actual forces encoun- 
tered at the tool tip. Alternatively, direct measurement of the forces at the end effector 
through sensors located on the jaws will provide accurate feedback of the exact forces 
to the surgeon. In addition, direct force measurement does not require friction model- 
ing in the system unlike indirect measurement techniques since the interaction forces 
in the direct method are obtained at the site of interaction. This is crucial since sur- 
geons typically palpate tissues to characterize them as healthy or unhealthy [17]. 

In this paper we present our results on design, development, and testing of an 
automated laparoscopic grasper that can provide tri-directional force feedback. Our 
design addresses the challenges mentioned above and focuses on the direct measure- 
ment of the tool/tissue interaction forces. In addition to grasping tissues, the sensors 
also allow the grasper to be used as a probe for probing soft-tissues. Through prob- 
ing tissue or an organ surface, a surgeon can diagnose a localized area and then im- 
mediately grasp and manipulate the area without having to change tools or loose track 
of the diagnosed area. This paper will discuss: (1) design and development of the 
laparoscopic grasper, and (2) evaluation of the laparoscopic grasper in characterizing 
artificial tissues. 



40 



Gregory Tholey et al. 



2 Materials and Methods 

2.1 Design and Development of the Laparoscopic Grasper 
with 3-D Force Measurement Capability 



The prototype of a laparoscopic grasper with 3-D force measurement capability has 
incorporated the advantages of our previous designs with the addition of a direct 
force measurement and a more compact modular design [14]. The design has in- 
cluded a grasping jaw with force sensors that have the ability to measure the grasping 
forces in three directions, namely F x , F , and F z . This improvement will allow for a 
more accurate measurement and feedback of the grasping forces at the tool tip com- 
pared to other designs in the literature. Another key advantage of our design is the 
modularity of the end-effector to convert between a grasper, cutter, and a dissector. 

The prototype consists of a DC 
motor (manufactured by Maxon), 
cable-driven pulley system, and 
grasping jaw with sensors (See 
figure 1, part A). This prototype is 
the initial design and future ver- 
sions will encapsulate the elec- 
tronics and route the wiring 
through the hollow tool shaft. The 
use of cabling allows for near zero 
backlash and low friction in the 
mechanism. Also, the prototype 
has the DC motor placed inline 
with the grasper at the back end of 
the prototype for a more compact 
and aesthetic design. The mecha- 
nism starts with a steel cable (transmission cable) that connects the motor pulley to 
the driving pulley and routes the cabling by 90 degrees using two ball bearings (See 
figure 1, part B). The driving pulley uses two additional steel cables (one for each 
jaw) to open and close the end effector. These cables travel from the driving pulley, 
through the fine tensioner and the shaft of the tool, to the jaw pulleys, and then back 
to the driving pulley, bypassing the fine tensioner on the return path. At the jaw pul- 
leys, the cables are wound in opposite directions to create the opposing motion of the 
jaws. Also, using the same driving pulley for the cables will ensure that each jaw 
rotates by the same amount relative to the tool. 

The tensioning of the device, as shown in figure 1, part B, is accomplished with 
three independent mechanisms to account for each of the steel cables. The cable and 
motor system is mounted on a sliding plate for a coarse tensioning of the driving 
cables. This tensioning is performed by adjusting two screws located at the back of 
the prototype under the DC motor. The fine tensioning of the driving cables is 
achieved through the fine tensioner, which contains a separate mechanism for each 
cable. The cable enters the mechanism and passes over two bearings. The distance 




Fig. 1 . Laparoscopic grasper prototype. 
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between these two bearings can be adjusted by a screw to increase or decrease the 
tension in the cable. Therefore, each cable can achieve an adequate amount of tension 
independent of the other. The third tensioning mechanism is used for the transmission 
cable between the motor and driving pulley. The motor mount is attached to the slid- 
ing plate and was designed with the ability to slide relative to the sliding plate. Two 
screws located at the back of the prototype on each side of the motor are used to ad- 
just the motor mount and the tension in the transmission cable. 

The end effector of the grasper consists of two jaws, one equipped with sensors 
and one without sensors (See figure 1, part C). The size of the jaw with sensors is 
approximately 0.6 inches wide by 1.75 inches long by 0.5 inches high. The jaw with- 
out sensors has similar dimensions with a reduced height. Both jaws have a rectangu- 
lar gripping surface with a surface area of 0.5 square inches. The jaw with sensors 
consists of piezoresistive sensors in the X and Y directions and a thin-film force sen- 
sor for the Z-direction (See figure 1, part C). While we realize the current dimensions 
of the jaw with sensors on this prototype would be slightly oversized for laparoscopic 
procedures, future versions using smaller and fewer sensors, for a more appropriate 
size, will be developed. The basic assembly of the upper jaw consists of a thin, steel 
gripping surface (0.05 inches), thin-film force sensor, and an upper assembly, which 
is assembled using an adhesive cyanoacrylate (see figures 2 and 3). The upper assem- 
bly contains 4 piezoresistive sensors (see figure 4) for the X and Y directions, top 
case, bottom case and cover. The four sensors are necessary for the two directions 
because they only measure compressive forces. In the future version of this prototype 
(which is currently under development), we have only three force sensors, one for 
each direction. 

The force sensors are positioned in the bottom case with space in front of the force 
sensing surface of each the sensor (see figure 5). The top case of the jaw is then 
placed over the bottom case and has four protrusions that fit into the space in front of 
the sensors and applies a force in the associated direction when the jaw is grasping an 



Jaw Design - Exploded view 




Fig. 2. Jaw assembly. 





Fig. 3. Prototype jaw with 
sensors incorporated in the 
grasper jaws. 




Fig. 4. Piezoresistive 
force sensor used in the 
prototype for measuring 
the lateral and longitu- 
dinal forces in the 
grasper jaws. 
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object (see figure 5). This is achieved by designing the jaw as a floating body relative 
to the rest of the grasper. The top case of the jaw is mounted directly to the entire 
mechanism onto the jaw pulley. The bottom case and top cover are placed above and 
below the top case respectively and secured to each other but not the top case. This 
forms an enclosure for the top case that is secured by 4 small screws. Also, there is a 
small gap (0.010 inches) on all sides between this bottom case and top cover assem- 
bly and the top case. This creates the “floating assembly” that is able to move with 
the grasped object and transmit the force exerted by the object onto the sensors. 
Therefore, the jaw with the sensors is capable of measuring forces in three independ- 
ent directions while maintaining a minimized design as necessary in laparoscopic 
tools. 



Piezoresistive Sensors (F.) 




Piozorosistivo Sensors (F.) 




Fig. 5. Cross section view of the jaw showing the piezoresistive sensors. 



The force sensors used in this grasper consist of one thin film sensor (model 402, 
manufactured by Interlink Electronics) for measurement of the normal force on the 
jaw and 4 piezoresistive sensors (FSS low profile force sensor, manufactured by 
Honeywell) for measurement of the lateral and longitudinal forces on the jaw. The 
thin film sensor is a force sensing resistor where an applied load changes the resis- 
tance through the sensor. The circuit used in conjunction with this sensor is a voltage 
divider with an op amp (see figure 6(a)). The range of this sensor is up to 25 N with a 
measured resolution of less than 20mN, as tested experimentally. The piezoresistive 
force sensor uses the circuit shown in figure 6(b). The potentiometers in this circuit 
allow for adjustment of the range and gain of the sensor. The range of the piezoresis- 
tive force sensor is 0 to 1500 grams with sensitivity of 0.12 mV/gram. The sensitivity 
shift for 0° C to 50° C is 5.5% of the full scale. 

The system comprises of the laparoscopic grasper actuated by a DC motor 
(model A-max36, manufactured by Maxon). The motor includes an incremental en- 
coder with a resolution of 0.18 degrees. The motor control is achieved by using the 
dSpace DS1103 controller board (manufactured by dSPACE, GmBH). We have de- 
veloped a program using the dSpace interface that allows a user to input a desired 
position of the jaws while also measuring the forces from the jaw sensors. We im- 
plemented a PD controller to control the position of the jaws. The control law is given 
by: 



T = K P (q d -q) +K Md-q) 



(1) 
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Fig. 6. Piezoresistive and thin film force sensor circuit diagrams. 



where T is the motor torque, K p and K v are the proportional and integral gains, q d and 
q are the desired and actual positions, and q d and q are the desired and actual veloci- 
ties. The inertial and friction effects in the mechanism do not have to be compensated 
for in the controller because a high gain is used and the grasper utilizes a direct force 
measurement method at the jaws. 



3 Experiments 

3.1 Characterization of Artificial Tissue Samples 

Our first experiment with the laparoscopic grasper was to evaluate the forces detected 
by the tool when grasping tissue samples of varying stiffness. For this experiment, we 
developed several artificial tissue samples made up of Hydrogel material. The Hy- 
drogel is created using a combination of polyvinyl alcohol (PVA) and polyvinyl pyr- 
rolidone (PVP) in the consistency of 90% PVA and 10% PVP. This solution was then 
caste into molds and subjected to several freezing/thawing cycles. The Hydrogel 
samples were cylindrical in shape with a diameter of 0.875 inches and a thickness of 
0.3 inches. Six cycles were performed and after each cycle a sample was removed. 
The tissues were numbered from 1 to 6 (in the order they were removed) with sam- 
ple 1 being the softest tissue while sample 6 was the hardest tissue. 

For the experiment, we selected three (1,3, and 6) of the Hydrogel samples that 
had a significant variation in stiffness such that they would be easily differentiated by 
direct exploration with one’s fingers. Each tissue sample was fully grasped by the 
tool and the normal force on the jaw was recorded. The tissue was grasped at the 
center of the jaw for each different sample to ensure all tissue samples were grasped 
in the same manner. In our current prototype, the tissue samples must be grasped at 
the center or back portion of the jaw. If the samples are grasped with the tip of the 
jaw, an incorrect force measurement will result due to reduced sensitivity at the tip. 
We are working towards developing a robust prototype in the second version of this 
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Tissue Characterization 




Fig. 7. Characterization of artificial tissue 
samples using the grasper (Magnified view of 
the forces shown in the figure). 



Tissue Characterization - Hysteresis 
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Fig. 8. The hard tissue sample showing 
hysteresis effects. 



device. The final position of the jaws was held constant and each of the samples was 
of the same thickness. Therefore, the force is the only variable in the setup. In Fig- 
ure 7, we plotted the norm of all three forces in three independent directions as sensed 
by the force sensors in the grasper. The norm of the forces was used to determine 
whether the tissue sample was hard, medium, or soft. As shown by the results in 
figure 7, the grasper can distinguish between samples of different stiffness. The soft 
Hydrogel sample showed a peak force of 0.2 N while the medium Hydrogel sample 
showed a peak force of 0.55 N. The hard Hydrogel sample showed a 0.9 N peak 
force, which was larger than the soft and medium samples. We also observed from 
Figure 7 that as the tissue is grasped, the force vs. displacement profile for the hard 
tissue has a larger force value throughout the compression of the tissue compared to 
the other two samples. As the tissue samples become stiffer, the required force to 
achieve the same compression also increases. Therefore, the grasper’s capability of 
differentiating between tissues of varying stiffness has been shown. 

Additionally, we observed hysteresis during the grasping of the tissue samples. 
Figure 8 shows the full grasping task (closing and opening the jaw) with the forces 
detected on the jaw. The grasper would close the jaw, hold the tissue for approxi- 
mately 10 seconds and then open the jaws. As shown in figure 8, there is some hys- 
teresis (up to 0.1N) in the force measurement. This hysteresis was observed to be 
caused by drift in the force sensor during the hold phase of the grasping task. How- 
ever, the hysteresis is mostly located at the point where the jaw begins to open. 



3.2 Direct vs. Indirect Force Measurement 

Our second experiment with the grasper involved the comparison of direct force 
sensing and indirect force sensing technique. Direct force sensing involves using a 
force sensor at the exact location where the measurement is desired while indirect 
force sensing might involve placing the sensor away from the location of the meas- 
urement. One example of indirect force measurement involves placing a strain gage 
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on the handle of a laparoscopic tool to measure the force at the tool tip through an 
appropriate calibration. The direct method in this experiment was using the thin film 
force sensor located on the jaw of the grasper to detect the forces. The indirect force 
estimation method involved calibration of the motor torque relative to the end effec- 
tor force. This calibration involved removing the bottom jaw of grasper and placing a 
force sensor (manufactured by JR3) under the upper jaw. By increasing the current 
supplied to the motor, the force exerted by the jaw on the force sensor increased and 
thus the relationship between the motor torque and end effector force can be deter- 
mined. The calibration rule after several trials was determined to be: 

Indirect force measurement: F. = 0.71 * V m (2) 



where Fj is the jaw force and V m is a controller input proportional to the voltage sup- 



plied to the motor. 



We used Hydrogel samples again to 
conduct this experiment that consisted 
of a simple grasping task. The soft Hy- 
drogel sample (sample 1) was grasped 
while recording the sensor readings, the 
motor torque, and the position of the 
jaws. Then the two methods of obtain- 
ing the force (direct and indirect) were 
plotted and compared. As shown in 
figure 9, there was a significant differ- 
ence between the two methods. While 
the indirect method showed a linear 
Fig. 9. Direct vs. indirect force sensing force curve from the calibration rule, it 
using the laparoscopic grasper. significantly overestimated the tissue 

grasping force by an order of magnitude 
compared to the direct force sensing method. This difference can be explained by the 
natural compliance of the artificial tissue, which is captured by direct force sensing. 
However, both methods arrive at approximately the same value of 0.18 N at the final 
jaw position. This experiment validates the requirement for direct force sensing. 
However, it should also be noted that indirect force sensing is significantly governed 
by the calibration procedure and the accuracy of the calibration equipment involved. 
In our experiments, while we did several tests to arrive at the calibration law in Eq. 
(2), we still observed a significant difference in the tissue grasping experiment de- 
scribed above. 




3.3 Soft-Tissue Probing Using the Laparoscopic Grasper 

Our final experiment using the grasper focused on the probing ability of the tool and 
the accuracy of the forces measured. While the thin film sensor plays a dominant role 
in sensing the grasping forces in both direct and indirect force measurement tech- 
niques, our final experiment was to measure the accuracy of lateral and longitudinal 
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Load Cell 



Fig. 10. Experimental setup to measure 
the longitudinal probing force. 



Piezoresistive Sensor Force Measurement 




Fig. 11. Longitudinal force measurement using 
the grasper compared to the actual forces meas- 
ured by a load cell during a probing task. 



forces measured by the automated laparoscopic tool. Figure 10 shows the experimen- 
tal setup to measure the normal probing force on a rigid object attached to the force 
sensor on one side and probed by the laparoscopic grasper on the other side. We used 
a rigid object for this experiment because we wanted to confirm the accuracy of the 
force measurement capability of the laparoscopic grasper by making certain that the 
probing force exerted on the rigid object is directly transmitted to the force sensor on 
the other side. 

A load cell (MLP-10, manufactured by Transducer Techniques) was used in this 
experiment. The grasper was then pushed forward towards the rigid object to simulate 
a typical probing task. As shown in figure 1 1 , the force measured by the load cell and 
the piezoresistive force sensor in the grasper jaws is of similar magnitude. There was 
very little error between the actual probing force (load cell) and the measured force 
by the grasper piezoresistive force sensors with the configuration in which they were 
places in the jaw (see Figure 2 for the jaw assembly). A similar experiment was per- 
formed to measure the lateral forces. Therefore, the capability of the grasper to be 
used to probe soft tissues or organ surface is a natural extension of this experiment. 



4 Conclusion 

This paper discusses the design, development, and testing of an automated laparo- 
scopic grasper with force measurement capability in three directions. The mechanism 
has low friction, near zero backlash, and is backdriveable. It also has a modular de- 
sign for interchangeability between various end effectors (cutter, grasper, dissector). 
Experiments were conducted to characterize artificial tissue samples, to compare the 
difference between direct and indirect force sensing methods, and to evaluate the 
grasper’s ability to measure probing forces. The results show the capability of the 
grasper to distinguish between Hydrogel samples of varying stiffness. The compari- 
sons of direct vs. indirect force measurements has shown that indirect force meas- 
urement can oveestimate the grasping forces and hence lead to significant error in 
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estimating the grasping forces of soft tissues. However, it is important to note that the 
grasping forces through indirect techniques can arise from the accuracy of the calibra- 
tion of the experimental setup and the changing contact conditions. The probing task 
has shown high accuracy when comparing actual probe forces to those measured by 
the jaw’s sensors. 

Future work with this prototype will consist of transmitting the forces measured 
by the grasper to a haptic interface device, such as the PHANToM, for providing 
accurate force feedback to the surgeon while manipulating animal tissues. In addition 
to adding the haptic interface, we will also attach the prototype of the end effector of 
a robotic arm and perform telemanipulation and diagnostic experiments on animal 
tissues. Also, future versions of this prototype will involve packaged electronics, only 
three force sensors, and a smaller design for usability in laparoscopic procedures. 
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Abstract. Exophthalmia is characterized by a protrusion of the eyeball. The 
most frequent surgery consists in an osteotomy of the orbit walls to increase the 
orbital volume and to retrieve a normal eye position. Only a few clinical obser- 
vations have estimated the relationship between the eyeball backward dis- 
placement and the decompressed fat tissue volume. This paper presents a 
method to determine the relationship between the eyeball backward displace- 
ment and the osteotomy surface made by the surgeon, in order to improve ex- 
ophthalmia reduction planning. A poroelastic finite element model involving 
morphology, material properties of orbital components, and surgical gesture is 
proposed to perform this study on 12 patients. As a result, the osteotomy sur- 
face seems to have a non-linear influence on the backward displacement. More- 
over, the FE model permits to give a first estimation of an average law linking 
those two parameters. This law may be helpful in a surgical planning frame- 
work. 



1 Introduction 

Exophtalmia is an orbital pathology that affects the ocular muscles and/or the orbital 
fat tissues [1]. It is characterized by a forward displacement of the eye ball outside the 
orbit (Fig. 1 (a)). This displacement, called protrusion, may leads to aesthetical prob- 
lems and to physiological disorders such as a tension of the optic nerve (dangerous 
for the patient vision) and the ocular muscles and/or an abnormal cornea exposition to 
the light. Exophthalmia can be due to four causes [1]: a trauma, a tumor, an infection 
or a disthyroidy. In our works, we mainly focus on the disthyroidian exophthalmia 
which is the result of an endocrinal dysfunction, as the Basedow illness. This pathol- 
ogy mostly leads to a bilateral exophthalmia since it induces a volume increase of the 
ocular muscles and/or of the fat tissues. 
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The classical treatment of exophthalmia is characterized by two steps. The first one 
aims to stabilize the endocrinal activities (by radiotherapy or surgery). The second 
step consists in a surgical reduction of the protrusion. The most efficient surgical 
technique is a decompression of the orbit [2] [3], that is to say an osteotomy of a part 
of the orbital bone walls via an eyelid incision. It leads to an increase of the orbital 
volume and thus offers more space to the soft tissues, particularly into the sinuses. To 
improve the backward displacement of the eye ball, some surgeons push on it in order 
to evacuate more of the fat tissues in the sinuses (Fig. 1 (b)). The whole procedure is 
critical since it can produce perturbations in visual functions (e.g. transient diplopia) 
and may affect important structures such as the optic nerve. 




Fig. 1 . (a) At the left; exophthalmia can be light (top), moderated (middle) or excessive (bot- 
tom). (b) At the right; the decompression is performed in the sinuses. The fat tissues then occu- 
pied those new cavities 



Up to now, the prediction of the results of an exophthalmia reduction is based on 
clinical observations [4] that gave the following average law: for a 1 cm 3 soft tissues 
decompression, a backward displacement from 1 mm to 1.5 mm is expected. Besides, 
few works dealing with the biomechanical modeling of the orbital soft tissues have 
been presented in the literature. [5] and [6] have developed a mass-spring model of 
the ocular muscles and the eye ball while [7] have proposed an elastic Finite Element 
(FE) model of the whole orbital soft tissues. Nevertheless, none of these models are 
applied to the exophthalmia reduction. 

In order to improve the exophthalmia reduction planning, a first study has been 
made by the authors [8] to predict the results of the orbital decompression. In this 
work, two complementary models have been presented to simulate the surgery. An 
analytical model that gave a good estimation of the volume decompressed in term of 
the eye ball backward displacement and a FE model that allowed to simulate the 
osteotomy and to compute the resulting backward displacement. In this paper, the FE 
mesh is used to study the influence of the osteotomy surface on the eye ball backward 
displacement. 
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2 Material and Methods 

In the previous work, the FE model was defined to correspond to a specific patient in 
order to estimate the variations between the simulation results and the clinical meas- 
urements, estimated on a pre-operative and a post-operative CT scans. To be as close 
as possible to the clinical set up for this patient, the morphology, the boundary condi- 
tions and the FE parameters were chosen to fit the data measured on this patient and 
taken as a reference. 

The first step was consequently to build the finite element mesh corresponding to 
this patient. To do this, the surface of the orbit was segmented to get the orbital bones 
surrounding the soft tissues. This segmentation has been done manually through the 
definition of a spline on each CT slice since those bones are thin and difficult to auto- 
matically separate from the rest of the tissues. The spline set is then extrapolated to 
give the three-dimensional surface of the bones. The geometry ant the volume of the 
muscles and the nerve can also be extrapolated with this process. 

From this 3D bone surface, the FE mesh can be generated. Since this geometry is 
complex, no software can be used to automatically mesh it. This process has been 
done manually using three-dimensional 20-node hexahedrons (quadratic elements). 
The resulting mesh (Fig. 2 (a) and (b)) is composed of 6948 nodes and 1375 ele- 
ments. 




Fig. 2. Manual generation of the orbital mesh, (a) left: in the patient skull and (b) right: with a 
dashed representation of the eyeball 

In order to simplify this model, the eyeball was considered as a rigid entity and was 
thus neglected. Since fat tissues occupy approximately 3/4 of the orbital soft tissue 
volume and since the fat tissues are predominant in the flow through the osteotomy, a 
homogenized material modeling the fat tissues has been chosen to simulate the orbital 
content. Clinical observations [1] describe the orbital fat tissues during a disthyroidy 
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as the combination of an elastic phase composed of fat fibers (mainly collagen) and a 
fluid phase composed of fat nodules saturated by physiological fluid. Consequently, a 
poroelastic material [9] has been used to model the intra-orbital soft tissues, using the 
FE software MARC© (MSC Software Inc.). 

The finite element properties of the model has been inspired by the literature ([10], 
[11] for the Young modulus and [12] for the Poisson ratio) and then adapted to fit the 
data measured on the pre-operative and post-operative CT scans of the patient studied 
in the previous work [8], This adaptation process led to a Young modulus value of 
20kPa and a Poisson ratio of 0.1. The two main poroelastic parameters, i.e. the poros- 
ity and the permeability, control the fluid behavior through the elastic phase and re- 
spectively take into account the fluid retention and the fluid pressure variation. Con- 
sidering the work of [13] for the inter vertebral disk, the orbital fat tissues permeabil- 
ity was set to 300mm 4 /N.s and the porosity to 0.4. 

Boundary conditions (Fig. 3) have been introduced to simulate the surgery: 

The constraint of the bone walls surrounding the soft tissues is translated by a nil 
displacement and a total sealing effect at the surface nodes. 



Fig. 3. Boundary conditions applied to the orbital mesh to simulate the decompression surgery 

Since the periost (the membrane around the soft tissues) remains after the wall os- 
teotomy, it is still constraining the fat tissue elastic phase. Nevertheless, the fluid 
phase is able to flow through this opening. Consequently, the osteotomy surface 
nodes are fixed in displacement while they are released in term of sealing effects. To 
study the osteotomy surface influence on this patient, four different osteotomies have 
been designed, in accordance with clinical feasibility. 

As the soft tissue volume increases, the pressure of these tissues increases too. This 
overpressure has been estimated at lOkPa by [14], This value was therefore applied as 
an initial pore pressure to all the FE mesh nodes. A relaxation time of 2s was then 
needed to reach pressure equilibrium. During this time, nodes located at the soft tis- 
sue/globe interface were free to move. 
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To simulate the force exerted on the eyeball by the surgeon, an imposed axial load 
has been applied to all the nodes located at the soft tissue/globe interface (via an ex- 
ternal controlled node to simply control this constraint). This imposed load was esti- 
mated clinically with a stiffness homemade sensor and gave an average maximum 
value of 12N. The load constraint is applied in 2s according to a linear ramp from 0 to 
12N. Since the eyeball is not modeled no contact analyze is required. 

The imposed load is maintained 3s to reach an equilibrium in the tissues and then it 
is released. To simulate the fact that the fat tissues stay in the sinuses after the de- 
compression, the osteotomy surface was set impermeable to avoid the fluid flowing 
back into the orbit. 

During this simulation, the eyeball backward displacement is estimated through the 
computation of the displacement of the external node where the imposed load was 
applied. 

This first work gave results close to the data measurements made for the patient 
studied particularly in terms of backward displacement estimation. Moreover, it ap- 
peared that the surface of the osteotomy had a preponderant influence on the eyeball 
displacement: not surprisingly, a larger osteotomy led to a greater backward dis- 
placement. 

Nevertheless, since this previous study has been made on only one patient, those 
results have to be confirmed. This is the aim of this paper which proposes a series of 
simulations on 12 patients with four different osteotomy surfaces. Fig. 4 shows those 
osteotomies. Their surfaces have a value of 0.8cm 2 , 1.7cm 2 , 3.4cm 2 and 5.9cm 2 . 




Fig. 4. The four osteotomies simulated to study the influence of the surface on the eyeball 
backward displacement 
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Since the first patient mesh already exists, an automatic meshing method has been 
used to generate the 1 1 new patients meshes (Fig. 5). This method is called the Mesh- 
Matching algorithm and is described in [15]. This algorithm is based on an optimiza- 
tion method that computes the elastic transformation between a reference mesh and a 
set of surface points. In this case, the reference mesh is the patient mesh that already 
exists while the set of surface points is given by the orbit segmentation for the 11 
other patients. By using this algorithm, the 1 1 new FE patient meshes are automati- 
cally generate in few minutes enabling to perform the FE simulations of the orbital 
decompression. The geometry of these meshes is various since the morphology of 
these patients are different. For example, the volume of the orbital soft tissues ranges 
from 18.1cm 3 to 31.4cm 3 . 




„ " 'Mesh-Matching Patient 2 
-Transformation 




Reference patient 




Patient 1 



Fig. 5. Transformation with the Mesh-Matching of the reference patient mesh in order to auto- 
matically generate the new patient meshes and to fit their morphologies. Here, two new patients 
with various morphologies are generated 



Since no post-operative CT scan was performed on these new patients, there was no 
possibility to analyze the behavior of the soft tissues during and after the surgery. 
Consequently, the FE parameters can not be evaluated for each of these patients. 
Though a variation of these parameters may exist clinically, the values used for the 
first patient were kept for these new patient models. 



3 Results 

Each of these simulations take roughly lh on a PC equipped with a 1.7GHz processor 
and IGo of memory. To estimate the influence of the osteotomy surface on the eye- 
ball backward displacement, a graph showing the relationship between those two 
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variables has been computed (Fig. 6). This graph clearly shows the non linear influ- 
ence of the surface on the backward displacement. A consequent increase of the os- 
teotomy surface is needed to get a moderate increase of the backward displacement. 
For example, an osteotomy surface 2 times larger leads to an increase of the eye ball 
backward displacement of 40%. 

Moreover, this graph points out the influence of the patient morphology (the vol- 
ume, radius and length of the orbit) on the backward displacement. Indeed, the curves 
presented in Fig. 6 are roughly ordered by the orbital volume, with the smallest vol- 
ume at the bottom of the graph and the largest orbital volume at the top. A difference 
of about 50% can be measured between the two extreme patients. Nevertheless, no 
generic law taking into account the patient morphology influence has been found to 
match the results of these simulations. 

Consequently, an average law has been computed to estimate the relationship be- 
tween the osteotomy surface and the eyeball backward displacement. This law is the 
equation of the average curve computed from the 12 patient curves. It is not a precise 
estimation of the clinical behaviour of the orbital soft tissues during the surgery since 
it is far from the extreme values as shown on the Fig. 6. This equation gives the back- 
ward displacement disp in function of the osteotomy surface .surf: 

disp = 1 . 1 * \n(surf) + 1.9 . ( 1 ) 



The following equation seems more useful to a surgeon that the 1cm 3 versus 
1/1. 5mm relation proposed by [4] . Indeed, it gives an estimation of the osteotomy 

surface, surf, to do in order to obtain a suited eyeball backward displacement, disp: 

disp— 1.9 



surf = e u 



( 2 ) 
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Fig. 6. Influence of the osteotomy surface on the eyeball backward displacement for the 12 
patients and average curve 
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4 Discussion 

The FE model presented in this paper seems to be able to answer the main question 
raised by the surgeons during the exophthalmia decompression: what osteotomy 
surface has to be done to get a suited backward displacement? The answer is esti- 
mated with equation (2). Indeed, the previous study has shown that the results given 
by the poroelastic model were close to the clinical measurements. Even if no meas- 
urements have been done for the 1 1 new patients, because of the absence of post- 
operative CT scans, it can be assumed that the results plotted in the Fig. 6 are not too 
far from the real soft tissues behavior. 

With this FE model, the non linear influence of the osteotomy surface seems non 
neglectable. Indeed, an important increase of this surface must be done to get an ef- 
fect on the eyeball backward displacement. Given that this surface cannot be ex- 
tended infinitely, since the orbital cavity is physiologically small, it could be interest- 
ing for a surgeon to know the relationship between the osteotomy surface and the 
backward displacement during the surgery planning step. The equation (2) gives that 
relationship in real time and therefore could help a surgeon to take into account the 
influence of the surface to choose and, if it is necessary, to increase the osteotomy 
size knowing its relative influence on the eyeball displacement. 

Nevertheless, it has been shown that equations (1) and (2) do not take in considera- 
tion the influence of the patient morphology. Those equations are consequently im- 
precise. To be more precise, this study has to be done on more patients and with 
comparisons with post-operative CT scans, to check the validity of the model results. 
Those new simulations would permit to take into account the influence of the orbital 
patient morphology and consequently to improve the estimation of the two equations 
presented above. 

Another important assumption has been made during the definition of the patient 
models. Indeed, the same FE parameters for the orbital soft tissues have been used for 
each patient. Clinically this assumption does not seem to be valid since the orbital 
soft tissues tension and the degree of the muscle fibrosis are heterogeneous and may 
vary from one patient to another. Here again, a most important study, with post- 
operative scans, would be able to take into account this fact and to determine the 
inter-patient FE parameter variations. 



5 Conclusion 

This paper has proposed an estimation of the influence of the osteotomy surface on 
the eyeball backward displacement in the context of orbital decompression. This 
evaluation is based on a poroelastic finite element model of the orbital soft tissues, 
presented in a previous study. It seems to point out that the relationship between the 
surface and the backward displacement is non linear (an important increase of the 
surface leads to a moderate increase of the displacement) and that the patient mor- 
phology may also influence the eyeball displacement. The equation (2) gives a first 
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real time estimation of the law linking the backward displacement to the osteotomy 
surface. This equation may be helpful for a surgery planning and may replace the 
actual law stating that for a 1 cm 3 soft tissues decompression, a backward displace- 
ment from 1 mm to 1.5 mm is expected. Nevertheless, the study will have to be per- 
formed on a more important patient set to verify and to improve this estimation. The 
Mesh-Matching algorithm will consequently be useful in the future since it gives an 
easy process to generate automatically new patient meshes. 

Future works will thus focus on a further study on a more important set of patients 
with post-operative scans to (i) try to determine the influence of the patient morphol- 
ogy on the backward displacement (and improve equation (2)), (ii) estimate the varia- 
tions of the FE parameters amongst different patients and (iii) validate the results of 
this study. 

As a perspective, the integration of the eyeball and the muscles in the FE mesh will 
be studied to take into account their actions on the soft tissues behavior. A contact 
fluid/structure model may consequently be developed and compared to the poroelas- 
tic model. 
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Abstract. Realistic tissue models require accurate representations of the proper- 
ties of in vivo tissue. This study examines the potential for tactile imaging to 
measure tissue properties and geometric information about subsurface anatomi- 
cal features such as large blood vessels. Realistic finite element models of a 
hollow vessel in a homogenous parenchyma are constructed in order to estab- 
lish a relationship between tissue parameters and tactile imaging data. A linear 
algorithm is developed to relate the tactile data to linearized tissue parameters. 
The estimation algorithm shows low errors in estimating the model parameters. 
A preliminary study on two porcine livers results in errors on the order of 20% 
in estimating the liver geometry. This result is promising given the small sam- 
ple size and parameter recording limitations of this preliminary study. Further 
work will reduce these sources of error and lead to in vivo testing with a mini- 
mally invasive tactile imaging scanhead. 



1 Introduction 

Realistic models for surgical simulation require accurate representations of the prop- 
erties of in vivo tissue. Recent measurements of organs such as the liver have used 
specialized apparatus to characterize mechanical properties adjacent to the organ 
surface [1,2]. This study examines the potential for tactile imaging to measure tissue 
properties and geometric information about subsurface anatomical features such as 
large blood vessels. In addition to informing surgical simulation, this technique may 
provide a fast, simple, and noninvasive means of intraoperatively locating vessels and 
diagnosing diseases such as cirrhosis that are characterized by changes in mechanical 
properties [3]. 

Tactile Imaging uses an array of pressure sensors to map the surface pressures that 
result from indenting the tactile imager into the surface of a soft material [4], This 
medical imaging modality quantifies the qualitative information provided by the hu- 
man sense of touch through palpation. Tactile imaging evolved as a means of detect- 
ing and characterizing pathologies in the human breast that manifest as areas of in- 
creased stiffness [5,6] and has progressed to include similar pathologies in organs 
such as the prostate [7,8]. Previous work [9,10] has resulted in algorithms for estimat- 
ing the parameters of stiff inclusions embedded in soft tissue, including lump diame- 
ter, depth, and modulus, as well as the modulus of the surrounding soft tissue. In vitro 
mean absolute errors (MAE) for these algorithms are 5-16%, while clinical assess- 
ment shows 13% MAE for lump diameter. 
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These algorithms can be readily extended to the problem of estimation of liver 
properties. Liver anatomy shows a macroscopically homogenous parenchyma punctu- 
ated by large branches of the hepatic vein (Figure 1). These vessels leave a signature 
of lower pressure on the surface pressure data collected in tactile imaging, and this 
information can be used to estimate parameters of the vessels and the tissue in which 
they are embedded. Following the approach of previous work on tactile signal inter- 
pretation [4,10], we first characterize the forward relationship between tissue parame- 
ters and the tactile signal using mechanical models (Figure 2). We then develop a 
linear algorithm for inverting the relationship to estimate tissue parameters from tac- 
tile images (Figure 3). The algorithm is then tested on porcine livers. 




Fig. 1. Cross-section through a porcine liver, showing large hepatic veins (indicated in white). 




Tactile Imager 
Scanhead 

Liver 

Parenchyma 

Vessel 



Fig. 2. Model of tactile imager and liver, with a homogenous parenchyma and a single round 
vein. 



2 Algorithm Development 

A closed-form relationship between the parameters of interest and the tactile signal is 
not feasible due to the large strains and contact interactions of the imaging process, so 
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Fig. 3. Imaging and signal interpretation approach. 



finite element method (FEM) models are used to generate the tactile images for 
ranges of the tissue parameters of interest. A typical cross-section through a liver is 
shown in figure 1. This main features of this cross-section are captured by the two- 
dimensional model shown in figure 2; a plane strain model is used to minimize com- 
putational efforts; the plane approximates the situation along the centerline of the 
cylindrical tactile scanhead. 



2.1 Mechanical Modeling 

We constructed finite element models based on figure 2, modeling the tactile scan- 
head on the device used in the experimental validation, a 16 x 16 array of capacitive 
sensors spaced 2 mm on center mounted on a section of a cylinder with a 38 mm 
radius. The location and orientation of the scanhead in 6 DOF was determined by a 
magnetic tracker (miniBird, Ascension Technologies, Burlington VT). The model 
represented the tissue as incompressible, isotropic, and linearly elastic. The interac- 
tion between the scanhead and the liver was assumed well-lubricated (i.e. frictionless) 
and the bottom of the liver was fixed to the substrate. 

The tissue parameters of interest are the background modulus, B, the tissue thick- 
ness, f, the vessel diameter, d, the vessel depth z, and the vessel pressure, V. The range 
for the modulus B of the tissue encompassed average moduli ranging from human to 
porcine livers [11]. The ranges for the geometry parameters were based on expected 
human anatomy [12]. The values of the parameters used in the thirty-six models cre- 
ated are specified in table 1 . The vessels were set in the middle of the tissue so that in 
the models, z = (t - d)/2. Preliminary studies indicated that the vessel pressure had 
insignificant effect on the tactile data for physiological ranges, and so the pressure 
was set to zero for the models used here. Tactile data was calculated for every 2.5 mm 
displacement of the scanhead for 40 mm to either side of the vessel, approximating 
experimental data collection. 

Table 1. Parameters for the finite element models constructed. 

Parameter Value 



B 

T 

d 



10, 12.5, 15 kPa 
40, 50, 60 mm 
5, 6.5, 8, 10 mm 
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2.2 Inversion Algorithm 

For each model (i.e. combination of tissue parameters), the calculated tactile pressure 
data at each 2.5 mm displacement was concatenated into a single row vector P. Simi- 
larly, the tissue parameters were assembled into a column vector G = [B t d z] T . We 
then characterize the problem as a linear inversion, and seek the transformation matrix 
A that minimizes the error e in G - AP + e. The relationship between the pressure in 
P and the parameters [B t d z], however, is not directly linear. Rather use the tissue 
parameters directly in G, we therefore look for functions of the parameters from 
which the parameters [ B , t, d, ~] can be calculated, but which are more linearly related 
to the tactile information. 

For example, if we approximate the tissue far from the vessel as a linear spring, the 
scanhead force is related to the equivalent spring constant Bit. This suggests that B 
and the function 1 It are approximately linearly related to the surface pressure data. 
Following similar arguments for the other parameters, we then use the input parame- 
ters G = [B l/td/z l/z] T . 



2.3 Finite Element Model Parameter Estimation Results 

For each of the model data sets, the other 35 models were used to generate the trans- 
formation matrix used in the estimation algorithm using the pseudo inverse A = 
G( P l P)' P 1 . The results of estimating the model tissue and vessel parameters are 
summarized in table 2. 

Table 2. Mean Absolute Error in estimating the underlying parameters of liver finite element 
models. 



Parameter 


Mean Absolute 
Error in Estimation 


Liver Modulus B 


0.8% 


Liver Thickness t 


9.3% 


Vessel Diameter d 


8.1% 


Vessel Depth z 


7.2% 



3 Preliminary Physical Testing 

3.1 Physical Data Collection 

The inversion algorithm was validated on two porcine livers from healthy 40 kg pigs 
harvested within one hour of sacrifice. The livers were immediately flushed with 
Heparin to minimize clotting, and perfusion with physiological saline solution at 36°C 
commenced approximately one hour later. Resulting tactile images are shown in fig- 
ure 4. 

Eight vessels were found in the two porcine livers, and multiple tactile images 
were made of several of the vessels, resulting in 14 usable maps for testing the inver- 
sion algorithm. For each set of tactile image data, tactile frames were collected every 
2.0 mm for a 40 mm linear region centered on the vessel. It was noted that the thin 
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Fig. 4. Tactile maps (averages of the spatially registered image sequences) of sections of por- 
cine liver lobe, showing decreased pressure over vessels. The left image shows one vessel 
spanning the width of the image, indicating the presence of a large vessel beneath the surface. 
The image at right shows two vessels running from left to right, with the upper one leaving a 
smaller impression in the tactile image (due to smaller size or greater depth). The images 
shown are approximately 80 mm x 40 mm. 



porcine liver lobes vary considerably in thickness, and even over a 40 mm section the 
thickness varied up to 15 mm. Since the images were obtained using a hand-held 
sensor the pressure applied varied for each frame. The average pressure across all the 
frames of interest was 18.2 Pa with a standard deviation of 10.3 Pa. 

For each set of data in turn, the transformation matrix was found using the other 13 
sets of data and the parameter estimation tested on the set in question. The actual 
parameters were recorded after tactile imaging was complete by dissecting the liver 
lobes and measuring the vessel diameter, depth from surface, and total tissue thick- 
ness. 



3.2 Porcine Liver Parameter Estimation Results 

Due to the small sample size, the background modulus of the livers studied was as- 
sumed to be the same across all samples. The results of estimating the geometry pa- 
rameters using our inversion algorithm are summarized in table 3. 

Table 3. Results of estimating the underlying parameters of porcine livers with large embedded 
veins. 



Parameter Mean Absolute 

Error in Estimation 



Liver Thickness t 
Vessel Diameter d 
Vessel Depth z 



20 . 0 % 

25.6% 

13.6% 
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4 Discussion 

The algorithm tests on the finite element models showed excellent accuracy, with all 
mean errors under 10% across the range tested. Estimation results from physical livers 
resulted in errors approximately twice those of the finite element models. This is not 
surprising, given the unconstrained data collection in the laboratory setup and the 
nonlinearities inherent in the tissue properties that are not captured by a linear algo- 
rithm. A key difference between finite element analysis and physical data collection is 
the large range of input pressures observed during physical data collection. This vari- 
able is controlled to better than 1% in the finite element analysis, but was observed to 
vary more than 50% in the physical data collection, due to human operation of the 
tactile imaging system. 

This variable input pressure range affects the data in two ways. First, due to the 
nonlinearity of the tissue modulus, we inadvertently probe the tissue at different effec- 
tive moduli. This precluded simple normalization of the frame information by the 
difference in the total applied force. This change in the effective background modulus 
will result in an incorrect measure of the tissue thickness, as we had assumed a con- 
stant background modulus. The wide range of input pressures also affects the way the 
tissue is probed in that as the tactile imager is indented further into the tissue, the 
surface pressure effectively senses deeper tissues [13], which results in information 
that our simple parameter system does not characterize. This problem of a wide range 
of applied forces can be alleviated by implementing bounds on the tactile data, using 
only data in a fixed range of total force, and signaling the user when they are operat- 
ing in the acceptable range. 

The preliminary experimental evaluation did not include estimation of the back- 
ground modulus because only two livers were available for testing. We note that pre- 
vious studies that estimated tissue modulus for in vitro models were successful, with 
5.4% MAE [10]. Because modulus was not estimated in the present study, the livers 
were assumed to have a constant background modulus. This assumption is reasonable 
as the subject animals were both healthy, the same age, and raised together. Any de- 
viations of the actual modulus will result in errors similar to those mentioned above 
for errors caused by a changing input force. 

With only eight vessels found in the two livers studied, the parameter range was 
sampled only sparsely. Obtaining extra maps on some vessels allowed for estimation 
of all map parameters, since we were not extrapolating for any one variable. Maps 
taken on the same vessel were not identical, and so represented different maps taken 
on vessels with similar parameters. Since these double maps were not identical, how- 
ever, they may have negatively affected the parameter estimation, by providing a 
different pressure signature for the same parameters. The duplicate maps were well 
spaced over the tissue thickness and vessel diameter, but were biased towards larger 
diameter vessels since the smaller vessels were difficult to image repeatedly, most 
likely due to temporary collapse. Therefore, although the total parameter range was 
spanned, the estimation of the vessel diameter was most likely adversely affected 
since the range spanned by the majority of the data was narrow, with more than one 
pressure profile representing the same diameters. 

The actual liver parameters were recorded after all maps were taken, by cutting the 
lobe perpendicular to the vessel along the line of data recording. Recording the pa- 
rameters this way is the most direct and readily available method, although it may 
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have contributed to inaccuracies in vessel parameter information. Since the cutting 
and data recording were done by hand, the planes of tactile data and dissection may 
be offset by a few millimeters. In this range, the vessel diameter and tissue thickness 
may vary as well. The vessel diameter may vary by up to a millimeter and the tissue 
thickness by twice that. The liver parenchyma also was prone to swelling in the cut 
plane. This is due to the natural tension that is present in the liver, maintained partly 
by the perfusion under which the data was recorded. Perfusion was necessary, how- 
ever, in order to maintain mechanical viability of the liver, so that despite the above 
sources of error, subsequent maps recorded on the same vessel record approximately 
the same conditions. These sources of errors largely affect the input parameters, and 
may adversely affect the apparent estimation by presenting incorrect information for 
the creation of the transformation matrix. These errors can contribute a relatively 
large error to the parameters in question, and so the above work should only be con- 
sidered a proof of concept for the use of tactile scanning to record liver vessel 
parameters. Within these constraints, the algorithm performed remarkably well in 
estimating the underlying parameters. 

Further experimental work should be conducted on livers after the main causes of 
error noted above are addressed. A method for regulating the force to a near-constant 
level should be implemented. This can be as simple as generating a specific sound for 
when the data collected is in a narrow range around the ideal input pressure [9], An 
improved method for measuring the geometry parameters without cutting the tissue 
should be employed. For further ex vivo tests using an accurate method such as MR1 
or CT scanning should be employed. A method to measure the parenchymal stiffness 
independently should also be utilized, so that the ability of tactile imaging as a means 
of recoding the background stiffness can be assessed. Ideally, human ex vivo livers 
will be tested before moving ahead to an in vivo setting using a smaller tactile imager 
in minimally invasive data recording. 
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Abstract. Mathematically describing the mechanical behavior of soft tissues 
under large deformations is of paramount interest to the medical simulation 
community. Most of the data available in the literature apply small strains 
(<10%) to the tissue of interest to assume a linearly elastic behavior. This paper 
applies a nonlinear hyperelastic 8-chain network constitutive law to model soft 
tissues undergoing large indentations. The model requires 2 material parameters 
(initial modulus, locking stretch) to reflect the underlying physics of deforma- 
tion over a wide range of stretches. A finite element model of soft tissue inden- 
tation was developed and validated employing this constitutive law. Ranges of 
the initial shear modulus and locking stretches were explored based on values 
found for breast tissue [17, 25]. Results of the model are shown with a lookup 
table containing third order polynomial coefficient fits. This work serves as an 
initial method to determine the unique material parameters of breast tissue from 
indentation experiments. 



1 Introduction 

Accurate mathematical descriptions of the mechanical behavior of soft tissues remain 
the limiting factor in the advancement of realistic medical simulations and non- 
invasive diagnostic tools. This is due to the complex nonlinear material properties of 
soft tissues when they undergo large mechanical deformations during minimally inva- 
sive procedures and diagnostic palpations. 

A phenomenological approach is implemented to realize the material parameters of 
soft tissues. Inverse finite element modeling (FEM) is used to fit mathematical ex- 
pressions in the form of a constitutive law to experimental data. Soft tissues are most 
often tested in an ex-vivo state with specimens of finite thickness under controlled 
loading and boundary conditions [8, 15, 18]. Selecting the appropriate constitutive 
law allows FEM to then be used to predict the tissue’s response to modes of deforma- 
tion not capable of being experimentally measured. We are interested in modeling the 
indentation of soft tissues by a rigid flat-ended cylindrical punch. 

This paper provides a method for determining an initial estimate of the material pa- 
rameters of soft tissue using the Arruda-Boyce constitutive model. A range of the two 
physically based material parameters of this model were explored based on data found 
in the literature on normal and pathologic breast tissue [17, 21, 22, 25]. The resulting 
force-nominal strain plots were fitted to 3 rd order polynomials whose coefficients and 
the resulting material parameters are presented. Others have attempted to model 
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breast tissue under uniaxial compression and assumed linear elasticity [3, 20]. Han 
assumed quasilinear viscoelasticity with an exponential elastic response and modeled 
the breast under plain strain conditions because he used a rectangular shaped probe on 
thick specimens [10]. 

Results indicated here should serve as a means for identifying an estimate of the 
physiologically based material parameters of the Arruda-Boyce model. 



2 Background 



There exists a well-defined analytic solution for indentation by a rigid flat punch 
assuming infinitesimal strains, frictionless contact, and a semi-infinite incompressible 
elastic half-space [14, 24], To account for the finite thickness in indentation experi- 
ments on cartilage, Hayes expanded the analytical solution to include a term that is 
dependant on the sample thickness: 



P = 



2aEw 



f 



a w'] 




(l) 



where P is the applied force, E is the elastic modulus, a is the radius of the indenter, w 
is the depth of indentation, and K is a dimensionless term to account for sample thick- 
ness (Fig. 1) [11]. Zheng created a finite element model using equation (1) to explore 
the effects of nonlinear geometry, namely large deformations up to 15% nominal 
strain, compressibility, and friction on the indentation of cartilage attached to a semi- 
infinite rigid half space [28]. Our goal is to further this approach by introducing both 
larger strains (-50%), and material nonlinearities into a FEM of soft tissue under 
indentation. 




Fig. 1 . Conceptual diagram of the soft tissue indentation model. 



It is well understood that soft tissues are viscoelastic, anisotropic, inhomogeneous, 
and have nonlinear force displacement characteristics [9]. To simplify the mathemati- 
cal analysis, many researchers assume an initial isotropy, local homogeneity, and 
study tissue deformations in the linear regime of <10% nominal strain [17, 19, 25]. 
However, typical surgical manipulations are often much larger than 10% nominal 
strain. It has been shown that at larger strains an elastic contrast exists between tissues 
of different pathologic states [17, 25, 27]. Therefore more accurate representations of 
soft tissue behavior are needed. 
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Holzapfel suggests that only biological materials and solid polymers (rubber-like 
materials) undergo finite strains relative to an equilibrium state [12]. Therefore it 
should not be surprising that the hyperelastic constitutive models developed for elas- 
tomers have frequently been used to study soft tissues [5, 8, 9, 13, 15, 16, 18, 23]. 
Hyperelastic materials are considered initially isotropic and exhibit a nonlinear instan- 
taneous response up to large strains. There are two ways in which the strain energy 
functions for hyperelastic materials are derived: one based on continuum mechanics 
and the other based on statistical mechanics. 

We have created a FEM with a hyperelastic constitutive model based on statistical 
mechanics. We describe that model and the predictions it makes for large strain in- 
dentations of pathologic breast tissues. 



3 Materials and Methods 

3.1 Creating and Validating the Finite Element Model 

Using commercial finite-element software (ABAQUS 6.3-1, HKS, Rhode Island), the 
present investigation created an axisymmetric rigid indenter model to analyze the 
indentation of soft tissue (Fig. 2). The model was validated against the analytical 
solution presented in equation 1 (with both K=1.0 and K>1.0) and compared to Zheng 
et al’s [28] finite element model under infinitesimal strains before adding a nonlinear 
constitutive law. 




Fig. 2. Axisymmetric FEM of rigid body indenter and soft tissue mesh with frictionless contact 
and sliding boundary conditions on the base. 

The indenter was modeled as an analytical rigid body with a flat-ended cylindrical 
shape 12 mm in diameter. Initially the tissue (cartilage) was modeled as a deformable 
meshed layer and was assumed to be linearly elastic, isotropic, and incompressible 
with Young’s modulus E=100 kPa and Poisson’s ratio v=0.4999 1 [28], Its mesh con- 
sisted of 4-noded hybrid quadrilateral axisymmetric elements (CAX4H), finely biased 
in the immediate regions underneath the indenter as shown in Figure. 2. The contact 
between the indenter and the tissue was modeled using a “contact pair” where the 
indenter was specified as “master” and the tissue as “slave.” The contact property was 
defined as frictionless so that the tissue could freely slip beneath the indenter. Known 



1 Due to ABAQUS’ limitations when v=0.5, an incompressibility of v=0.4999 was used. 
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displacements were then applied to the reference node of the indenter that was ini- 
tially oriented on the surface of the tissue. The reaction forces generated by the FE 
simulation were recorded and plotted against strain. 

To validate the model against the case of the true analytical solution of a rigid flat 
punch (equation 1 where K=1.0), the boundary condition of the bottom surface is free. 
The tissue was unconstrained in the lateral direction and an aspect ratio (indenter 
radius to sample thickness) of 0.1 was used to approximate the semi-infinite elastic 
half space. Results from the FEM under 0.1% nominal strain are within 3.3% of the 
analytical form of the solution. 

The FEM was modified to validate against Hayes’ analytical model (k> 1.0 in 
equation 1). Fixed boundary conditions simulated the attachment of cartilage to a 
rigid bony layer [11]. Two aspect ratios were tested (0.2 and 1.0) to a strain of 0.1%. 
Less than 2% error occurs when the FEM accounts for finite tissue thickness and is 
compared to both Hayes’ analytical solution and Zhang’s FEM results at nominal 
strains of 0.1% (Table 1). 

Table 1. A comparison of the kappa value between our model (Liu) and the analytical solution 
of Hayes and the FEM of Zheng for 2 different aspect ratios. 



Aspect ratio 


Model 


K 


% Error 


a/h = 0.2 


Liu 


1.260 


- 


Zhang 


1.244 


1.25 


Hayes 


1.281 


-1.67 


a/h = 1.0 


Liu 


3.564 


- 


Zhang 


3.590 


-0.72 


Hayes 


3.609 


-1.24 



After the model was validated assuming linear elasticity under infinitesimal strains, 
the model was changed to simulate soft tissue indentation tests containing both geo- 
metric and material property nonlinearities. Experiments on breast tissue indentation 
found in the literature allow for frictionless contact between both the indenter and the 
specimen, and between the specimen and the testing surface [17, 25], Thus, the 
boundary condition on the bottom surface of the tissue in the model was uncon- 
strained in the lateral direction. To compare to the experimental breast tissue data, the 
indenter size was changed to 4 mm in diameter and aspect ratios of 0.5, 1.0, and 1.5 
were created [17, 25]. The indenter fillet was increased to 2xl0" 4 mm to allow the 
tissue to be compressed to 50% nominal strain. The mesh bias was refined until a 
model of each aspect ratio could reach the set strain of 50%. Strain rate tests per- 
formed on breast tissue suggest that viscous effects can be neglected [17, 25]. Hypere- 
lastic nonlinear material parameters were added to the material property definitions of 
the tissue in the indentation model. Specifically, the Arruda-Boyce constitutive law 
was selected. 



3.2 The Arruda-Boyce Constitutive Model 
3.2.1 Motivation 

The continuum mechanics approach for developing hyperelastic strain energy func- 
tions are empirical, need more than one experiment to realize their many material 
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parameters, and have a limited strain region over which their results are applicable. 
Although higher order models fit the data well, they are complex, computationally 
expensive, and unstable at high stretches [1,4]. Despite these difficulties they are still 
widely used to describe the behavior of soft tissues [5, 8, 9, 13, 15, 16, 18, 23]. 

The statistical mechanics treatment of rubber elasticity (Langevin chain statistics) 
model the material chain segment between chemical cross-links as a rigid link with 
set length [4], The stress- strain behavior is governed by changes in configurational 
entropy [2]. The end result reflects the underlying physics of macroscopic deforma- 
tion from microscopic components. In particular the Arruda-Boyce model is an 8- 
chain network model where only two material parameters (the rubbery initial modulus 
and the limiting chain extensibility) are needed to describe the behavior of a material 
over a wide range of stretches given limited test data. This model lends itself ideally 
to that of soft tissues because the polymer chains mimic their main constituents: col- 
lagen and elastin fibers. 



3.2.2 Development 

A convenient form of the Arruda-Boyce strain energy function, U, is found by taking 
a series expansion of the inverse Langevin function to the 5 th order: 



U=ni-Ml[- 3 ')+ 4 

i=l4„ L> 






■ ln./„ 



( 2 ) 



where 



c =— c =J-c =JLc =-^c = -™- 
1 2 ’ 2 20 ’ 3 1050 ’ 4 7000 ’ 5 673750 



(3-7) 



fi = nK9 

2„,=V(V (8-10) 

/j =tr(A ) = 2f +2^ +2^. 

Here p is the initial rubbery shear modulus, n is the chain density, 0 is the tempera- 
ture, K is Boltzmann’s constant, 7 m is the limiting chain extensibility (locking stretch), 
N is the number of rigid links, 1 is the first deviatoric strain invariant, and L are the 
principle stretches. D is a temperature dependant material parameter related to the 
bulk modulus, and J el is the elastic volume ratio. For incompressible materials, J el = 1 
and the second term in equation 2 drops out. 

Due to symmetry each chain’s stretch is shared equally amongst all of the chains 
and an initially isotropic configuration can be assumed. We can therefore relate the 
microscopic chain length to the macroscopic principal stretches via: 

\hain = (aMT + A )=J- ( 1 1 ) 

where 2 chain is the chain stretch, l is the current chain length, and l 0 is the initial chain 
length. The locking stretch, A, m , is the value of the chain stretch when it reaches its 
fully extended state, and can be determined from a simple tension or compression 
experiment assuming incompressibility (A.jA.jA^l). Modeling a uniaxial compression 
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state and noting from the literature that breast tissue drastically increases its stress 
response at strains on the order of 30% nominal strain for normal tissue and 10% 
nominal strain for cancerous tissue, locking stretches were calculated to be 1.05 and 
1.01 respectively. 

Literature on breast tissue material property measurements of varying pathology 
suggest that initial elastic moduli for cancerous tissue is between 3 and 7 times that of 
normal tissue [17, 21, 22, 25]. For an incompressible material Possion’s ratio is 0.5 
and the elastic modulus is equal to 3 times the shear modulus. Given initial elastic 
moduli reported in the literature of 33 kPa and 100-186 kPa for normal and infiltrat- 
ing ductal carcinoma respectively at 5% nominal strain, we explored initial shear 
moduli between 1 kPa and 150 kPa. 



3.3 Applying the Nonlinear Constitutive Law to the FEM 

Using the proposed FE model, we chose to model eight different combinations of the 
two material parameters of the Arruda-Boyce constitutive law over three different 
aspect ratios (a/h = 0.5, 1.0, 1.5): p=l, 5, 10, 20 kPa with = 1.05, and p=30, 60, 
100, 150 kPa with 7„ m = 1.01. The models were deformed to 50% nominal strain and 
the displacement and reaction force in the axis of deformation were recorded. 

The Arruda-Boyce model is typically used for very large strains (i.e. tensile strains 
> 200%). We are only interested in compressive strains on the order of 50%. We can 
assume the effects from the higher order polynomial terms are therefore negligible. 
Third order polynomials were fit to the force-nominal strain curves generated from 
our FEM analysis. The coefficients of these polynomials can be compared to experi- 
mental data to estimate the material parameters of the substrate under study. 



4 Results 

The force-nominal strain responses of the FEM with different values of the initial 
shear modulus and locking stretch material parameters of the Arruda-Boyce constitu- 
tive model are shown below in Figure 4. Eight values of the initial modulus ranging 
from 1 kPa to 150 kPa were modeled with two values of the locking stretch (1.01 and 
1.05) based on breast tissue data found in the literature as previously stated. Three 
aspect ratios were modeled to account for different sample thickness and indenter 
geometry. Typical computation times for 50% strain on a Pentium III computer were 
on the order of 140 seconds. The coefficients of third order polynomials fit to the 
model’s response and their R 2 values are shown in Table 2. 



5 Conclusions and Future Work 

For realistic medical simulations to become a practical reality the acquisition of bio- 
mechanical information and efficient computation must be achieved. The latter is left 
to the many researchers working on deformable meshing techniques [6, 7, 26]. It was 
the intent of this paper to focus on uniquely characterizing the complex nonlinear 
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Table 2 . The material parameters of the Arruda-Boyce model and their resulting third order 
polynomial fit coefficients for the force-nominal strain responses of soft tissue indentation. 



a/h = 0.5 



R A 2 


A 


B 


C 


lambda 


mu (kPa) 


0.9997 


1.48 


0.53 


0.31 


1.05 


1 


0.9997 


7.46 


2.72 


1.56 


1.05 


5 


0.9997 


14.80 


5.28 


3.11 


1.05 


10 


0.9997 


29.88 


10.83 


6.25 


1.05 


20 


0.9997 


57.43 


20.83 


10.90 


1.01 


30 


0.9997 


113.81 


41.21 


21.69 


1.01 


60 


0.9997 


187.01 


67.21 


35.87 


1.01 


100 


0.9996 


275.91 


98.22 


53.38 


1.01 


150 



a/h = 1 .0 



R A 2 


A 


B 


C 


lambda 


mu (kPa) 


0.9995 


1.20 


0.43 


0.23 


1.05 


1 


0.9995 


5.92 


2.10 


1.16 


1.05 


5 


0.9995 


11.73 


4.16 


2.30 


1.05 


10 


0.9995 


23.16 


8.17 


4.58 


1.05 


20 


0.9995 


43.93 


15.43 


7.93 


1.01 


30 


0.9995 


86.04 


29.99 


15.71 


1.01 


60 


0.9995 


140.78 


48.72 


25.93 


1.01 


100 


0.9995 


207.65 


71.35 


38.55 


1.01 


150 



a/h = 1 .5 



R A 2 


A 


B 


C 


lambda 


mu (kPa) 


0.9995 


1.19 


0.44 


0.22 


1.05 


1 


0.9995 


5.85 


2.17 


1.07 


1.05 


5 


0.9995 


11.56 


4.27 


2.13 


1.05 


10 


0.9995 


22.75 


8.38 


4.22 


1.05 


20 


0.9994 


42.98 


15.75 


7.31 


1.01 


30 


0.9994 


83.87 


30.52 


14.44 


1.01 


60 


0.9994 


136.84 


49.49 


23.79 


1.01 


100 


0.9994 


201.47 


72.42 


35.31 


1.01 


150 



behavior of soft tissues with a simple mathematical model given limited experimental 
data. We implemented such a model in finite element simulations to predict the be- 
havior of soft tissues undergoing large indentation deformations across various testing 
geometries. The model was validated in the linear regime against analytical solutions 
and another FEM. The force-nominal strain results of the model can be used to esti- 
mate the material parameters of the soft tissue of interest. 

An axisymmetric finite element model with frictionless contact and boundary con- 
ditions was created employing the hyperelastic Arruda-Boyce constitutive model. 
Unlike similar constitutive laws formulated from continuum mechanics, this statistical 
mechanics based model was chosen because its two material parameters have a physi- 
cal interpretation that can be directly related to the constituent make-up of soft tissues 
(collagen and elastin fibers). 

Most of the soft tissue data published in the literature only apply nominal strains of 
up to 10%. At these low strains, the Arruda-Boyce model reduces to the linear elastic 
Neo-Hookean form and fits the data well. At strains where the usefulness of these 
elastic models is of limited value the Arruda-Boyce model continues to predict the 
nonlinear behavior of the soft tissues. 
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Fig. 3. Force versus nominal strain results for the FEM with a/h=0.5 (top), a/h=1.0 (middle), 
and a/h=1.5 (bottom) across various initial shear moduli and locking stretches (L m =1.05 (left) 
and (X m =1.01 (right)) up to 50% nominal strain. 



Wellman has collected some indentation data on pathologic breast samples at lar- 
ger strains (>35% nominal strain). A future study will analyze this data and compare 
it to nonlinear FEM simulations to determine the unique material properties of the 
tissue. The lookup tables presented in this paper will be used to obtain approximate 
initial values for the material parameters. An iterative process will then ensue where 
the models’ results will be compared to the large strain data via the employment a 
nonlinear search scheme minimizing the sum of squares error. With an educated esti- 
mate of the initial value of the material parameters, convergence of a unique set of 
material parameters can be quickly obtained to within the standard error of the data 
collected. 

A preliminary set of both normal glandular and cancerous data are plotted together 
with the corresponding FEM results in Figure 5. Fitting a third order polynomial to 
the data in Figure 5 suggests that an estimate for the normal glandular tissue material 
parameters are on the order of X m = 1 .05 and u= I kPa (a/h=1.0) and for infiltrating 
ductile carcinoma A. m = 1.01 and p=30 kPa ( a/h= 1.5). 
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Fig. 4. Preliminary large strain indentation data plotted for normal breast tissue with a/h=1.0 
(left) and infiltrating ductile carcinoma (right) against FEM results with a/h=1.5. 



It is clear from this preliminary work that the model needs further development. 
The tissue appears to have a lower locking stretch than the Arruda-Boyce model pre- 
dicts. This is most likely because the Arruda-Boyce model assumes an initial stress- 
free reference state, where as in real tissue this does not exist due to hydration and 
tension in the fibers. Accounting for this non-zero initial stress state should bring the 
model into closer agreement with the data and is currently being developed. 
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Abstract. A 3D biomechanical model of the tongue is presented here. Its goal 
is to evaluate the speech control model. This model was designed considering 
three constraints: speech movement speed, tongue movements and tongue soft 
tissue mechanical properties. A model of the tongue has thus been introduced 
taking into account the non linear biomechanical behavior of its soft tissues in a 
large deformation analysis. Preliminary results showed that the finite element 
model is able to simulate the main movements of the tongue during speech 
data. It seems that the model may be used in the estimation of glossectomy im- 
pacts on patient speak and on different surgery approaches. 



1 Introduction 

The A 3D biomechanical model of the tongue is developed to serve as a tool for fu- 
ture evaluations of speech motor control models. This model was designed consider- 
ing three constraints, crucial for speech tasks. Firstly, speech movements are rapid 
(around 50ms) and their time variations determine important cues for speech percep- 
tion: the model must therefore handle dynamics [1] and be able to produce fast ges- 
tures. Secondly, tongue movements are resulting from the interaction between muscle 
commands, muscle anatomy, tongue tissues mechanical properties and external fac- 
tors such as the teeth, palate or pharyngeal walls. An accurate representation of 
tongue muscular structure, a complete description of vocal tract and a realistic ac- 
count of tongue elasticity are therefore required to design the model. Tongue soft 
tissues feature non linear mechanical properties, mostly due to muscle fibers inter- 
weaving. Moreover, as reported in [2], the relative tongue deformations can reach a 
rate of 200% in compression. The model has then to handle non linear mechanical 
properties and to assume a large deformation framework. Model design and simula- 
tions were carried out with the use of the Finite Element package ANSYS™. 

A basic obstacle encountered in the design of such a model arises from the lack of 
experimental data giving quantitative information about the mechanical properties of 
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tongue tissues. This is the motivation of this paper, in which the tongue model will be 
first briefly described before the experimental and computational procedure will be 
presented that we have set up to quantitatively measure the elastic properties of a 
human tongue. 



2 Design of the Mesh 

The Finite Element Method [3] was used for the discretization of the equations of 
motion. This method requires the design of a mesh, based on biological data provid- 
ing an accurate description of the tongue. The Visible Human Project®’ s female data 
set, in parallel with morphological studies of tongue musculature, was used to create 
the original mesh (Figure 1). This process is described in details in [4] while the defi- 
nition of internal muscular structures is explained in [5]. 

In order to facilitate further quantitative comparisons of simulations with experi- 
mental data collected on real speakers, the external geometry of the initial mesh was 
matched to the 3D tongue shape of a human speaker tongue shape, while preserving 
its internal structure at rest. Reference data consisted in the association of MRI im- 
ages, CT scans (for bony structures like jaw and hyoid bone), and plaster cast of the 
palate optically scanned to digitize the information. 




Fig. 1. 3D Finite Element mesh (Front view, Side view. Sagittal view) 



3 Mechanical Modelling 

3.1 Simulation Framework and Choice of the Material 

Tongue models previously published by [6], [7], [8], [9], and [TO] were based on 
small deformation framework hypothesis. However, [2] reported that tongue defor- 
mation can reach a rate of 200% in compression and 160% in elongation during 
speech tasks, which is totally out of the small deformation framework. This is why a 
large deformations framework was assumed here. 

Biological soft tissues can feature non linear, anisotropic and viscoelastic me- 
chanical properties [11]. The non linearities are accounted in the model, using hyper 
elastic material. A material is said to be hyper elastic if we can find an energy func- 
tion W so that its derivative respect to the strain E equals the stress tensor S: 
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5 = 



dW 

dE 



E = -(F t F-I ) 
2 



( 1 ) 



E is the lagrangian deformation tensor, F is deformation gradient and I is the identity 
matrix. 

For our model, the energy function W was approximated by a 5 coefficients 
Mooney-Rivlin material, so that: 

W = a 10 (/| - 3) + a w (/ i - 3) 2 + « 01 (/ 2 - 3) + « 02 (/ 2 - 3) 2 + «„(/, - 3 )(/ 2 - 3) . (2) 



where II and 12 are the first and the second invariant of the deformation tensor E. In a 
first approximation, the model was assumed to be isotropic. 



3.2 Evaluation of Tongue Mechanical Properties 

An indentation experiment was provided on a fresh human cadaver tongue (Figure 2, 
left). The indenter is a 14mm radius cylinder, applying increasing loads locally on the 
tongue. The vertical displacement of the indenter is measured as a function of the 
applied load. 




Fig. 2. Experimental indenter (left) and its numerical model (right) 



Unfortunately, these data can’t provide direct results for the constitutive law of the 
material. In order to get to the stress/strain relationship (i.e. the Mooney-Rivlin coef- 
ficients) from them, the indentation experiment was numerically simulated, with a 
finite element model of the lingual tissues laid on a table and in contact with an in- 
denter (note that figure 2 only plots a quarter of the simulation: axisymetric assump- 
tion). A first series of 5 experiments was done in the front part of the tongue, and 
another series was done in the rear. Loads are gradually applied from 0 to IN, with 
0.1N steps. During relaxation, loads were applied from IN to ON with 0.2N steps 
(Figure 3). 
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Serie 1 




Fig. 3. Example of results for a set of 5 experiments (Load experiment) 

It should be noted that, due to tongue plasticity, tongue shape after relaxation is 
not the same as initial tongue shape. Thus, after each relaxation, the tongue was 
manually reshaped so that initial conditions can be supposed to be essentially the 
same for each experiment. The mucosa covering the surface of the tongue has me- 
chanical properties different from the muscles. Hence, two additional series of meas- 
urements were carried out after removing the mucosa. In summary, two series of 
experiments provide data on muscular tissues only whereas the two other series pro- 
vide data including the mucosa. 

Then, the determination of the constitutive law was done in two steps. A first 
model of the numerical indenter with a single layer structure was designed to com- 
pute the Mooney-Rivlin coefficients for the anterior and posterior parts of the tongue 
separately. In a second step and given the coefficients for muscle tissue computed 
previously, the mucosa constitutive law was computed, once for the front of the 
tongue and once for the rear. Then an iterative process was used to find the Mooney- 
Rivlin coefficients so that simulated displacements fit with experimental displace- 
ments. For each of the four sets of experiments, one constitutive law was computed 
and the Mooney-Rivlin coefficients were calculated using the mean of the 5 dis- 
placement/force measured. An iterative algorithm, based on a dichotomic conver- 
gence process, was designed to compute these coefficients. 



3.3 Boundary Conditions and Equation Solving 

The model is assumed to be incompressible. The Poisson ratio I should then be as 
close as possible to a 0.5 value. However, increasing the i value dramatically in- 
creases computation time. Consequently, different values of I were tested, ending up 
by choosing a 0.49 value for f. This value corresponds to a volume variation lower 
than 2% with reasonable computation times. The gravity value was fixed to 
g=9.81m.s-2 and the density to 1000kg/m3. The corresponding mass value was there- 
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fore around 140g. To take into account the visco-elasticity of tongue tissues, the 
Rayleigh damping model was assumed, and damping coefficients were chosen to 
reach critical damping. 

The forces are applied along macro-fibers defined by a series of nodes along the 
edges of the elements. Forces are mainly concentrated at the ends of fibers and aim at 
reducing the fiber length when the muscle is activated. However, the fibers are ini- 
tially curved and, in order to account for the fact that muscular activation tends to 
straighten fibers, forces depending on the fibers curvature are distributed along the 
fibers (cf. [10]). 

These equations are solved by the FE package ANSYSTM which computes the 
constraints, the deformations and displacements at each node of the mesh. 



4 Results 

Figure 4 plots the experimental results and the estimated constitutive law for the rear 
part of the tongue without mucosa. The dots are experimental measurements and the 
curve represents the constitutive law (al0=0.42kPa and a20=12.5kPa). The maximum 
gap between experimental and simulated data is 0.22mm. 







Fig. 4. Experimental results and estimated constitutive law for the rear of the tongue without 
mucosa 

Once the constitutive law of tongue soft tissues was calculated, it was implemented in 
the biomechanical tongue model. The model was then tested with force levels com- 
patible with speech movements. Here is an example of deformations obtained by 
specific muscle activation. The duration of each simulation is 120ms and the forces 
are applied as a step function for all simulations. Tongue shapes in the mid-sagittal 
plane, at rest and at the end of the simulation, are shown in Fig 5 for the simultaneous 
activations of the Posterior Genioglossus (6N), and the Transversalis muscle (2N). 
Expected results are a back to front movement of the tongue with an elevation of 
tongue dorsum toward the palate. Our simulations are in agreement with these predic- 
tions. 
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Fig. 5. Model in rest position (left) and final position (right) obtained with parallel activation of 
Genioglossus posterior (6N) and Transversalis (2N) 



5 Conclusion and Perspectives 

Preliminary studies on each muscle activation on tongue shape show that main trends 
of deformations are handled by this model within durations and with level of forces 
that are fully compatible with speech data. 

Impact of each muscle on tongue shape will be evaluated more precisely in future 
works. This model will be evaluated by comparison of our simulations with meas- 
urements of the shape of the vocal tract in the mediosagittal plan. The impacts of 
these muscle activations are compared to MRI measurements of tongue shapes for a 
selected number of sounds, for which EMG activations were published in the litera- 
ture. 

This biomechanical tongue model is to be used for pathological applications as a 
tool to help surgeons planning medical gestures during glossectomy or sleep apnea 
surgery. In both cases, this model could be used to estimate the impact of surgery 
gestures on patient’s speaking abilities, and if necessary help surgeons to plan a dif- 
ferent surgery strategy, less disturbing in terms of speech production abilities. 
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Abstract. The knee joint is a key structure of the human locomotor system. 
Any lesion or pathology compromising its mobility and stability alters its func- 
tion. As direct measurements of the contribution of each anatomical structure to 
the joint function are not viable, modelling techniques must be applied. The 
present study is aimed at comparing cruciate ligaments models of different 
complexity using accurate parameters from MRI and 3D-fluoroscopy of a sin- 
gle selected subject during step up-down motor task. The complexity of the 
model was not very relevant for the calculation of the strain range of the cruci- 
ate ligaments fibres. On the other hand, three-dimensionality and anatomical 
twist of the modelled fibres resulted to be fundamental for the geometrical 
strain distribution over the ligament section. 



Introduction 

The knee plays a fundamental role in determining the human locomotor ability. Any 
alteration of its anatomical structures can compromise its function. The development 
of effective methods for surgical reconstruction and rehabilitation is of great clinical 
interest, regarding both joint replacement and surgical reconstruction of the main 
anatomical structures. This interest is demonstrated by the 259000 total knee re- 
placements, 25000 ligaments reconstructions and 15000 other repairs of the knee 
performed in the USA in 1997 as reported by the American Association of Orthopae- 
dic Surgeons (AAOS). For the development of these procedures, an accurate knowl- 
edge of the mobility and stability of the whole articular structure, as well as of its 
different anatomical sub-units, is necessary. The need for this deeper knowledge led 
to a bulk of in-vitro and in-vivo studies, which allowed to clarify several aspects of 
the physiological behaviour of this complex joint. In-vitro testing allows to directly 
observe and measure different aspects of joint mechanics, but not in physiological 
conditions. During its normal function, the knee lets the shank move with respect to 
the thigh, maintaining the stability of the structure under articular load and torque. 
These are the result of several contributions: the inter-segmental contact load, liga- 
ment tensioning, loads applied by the muscles, the inertia of body segments. All these 
contributions are strongly dependent on the analysed motor task, as well as on the 
physical characteristics of the subject. Thus, if we want to quantify the contribution of 
each anatomical structure in determining the physiological function of the knee, mod- 
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elling is the only possible solution, as direct measurements cannot be performed. 
Further more, a knee model can be a useful teaching and clinical tool for the investi- 
gation of the function of the analysed joint. 

The problem of knee modelling has been approached at different levels of com- 
plexity. Two-dimensional models were designed in order to investigate the role of the 
cruciate ligaments in simple conditions, such as isometric quadriceps contraction 
[1,2]. Three-dimensional models, including articular surfaces and ligaments, were 
also proposed. Even these more complex models were applied in conditions far away 
from those of the physiological knee [3-6]. The natural evolution of this approach is 
inserting the model into a context which allows to evaluate the boundary conditions of 
the knee-structure during the performance of a simple task of daily living [7], Even if 
the model is designed properly for the application devised, its potentials can be nulli- 
fied by the effect of errors within the definition of subject parameters and during the 
acquisition of experimental inputs. In previous modelling attempts, these errors were 
due to discrepancies in the origin of parameters and inputs, which were often obtained 
from different and non-homogeneous subjects. 

In order to avoid this possible source of error, in this paper, different cruciate liga- 
ment models were compared using parameters from a single selected subject analysed 
as accurately as possible. The specific geometry of articular surfaces and ligaments 
insertions were reconstructed using the three-dimensional reconstruction of seg- 
mented bone and soft tissues, obtained from Magnetic Resonance Imaging (MRI). 
The specific accurate kinematics was obtained from cine-fluoroscopic images of a 
step up-down motor task. Cruciate ligaments models of different complexity: from the 
simple bi-dimensional untwisted one to the more realistic three-dimensional twisted 
with circular insertion were compared. The aim was to select the best compromise 
between accurate anatomical description and model simplicity for the investigation of 
knee biomechanics. 



Material and Methods 

Overview. A subject-specific model of the right knee of a young male living subject 
(height 168 cm, weight 62 kg, and age 30 years) was developed from a high resolu- 
tion MRI data set. Three-dimensional outer surfaces of the biological structures of 
interest were generated. 

The subject performed step up-down with the knee under analysis inside the 
fluoroscopic field of view. The accurate 3D pose of the bones was reconstructed by 
means of single-plane lateral 2D fluoroscopic projections and relevant models previ- 
ously obtained. 

The cruciate ligaments fibres were modelled with six geometrical equivalents and 
relative fibres strain compared: 2D, 3D with rectangular insertions, 3D with circular 
insertion, each twisted and untwisted. 

The MRI Data Set. A data set of high resolution MRI images was collected with a 
1.5T Gemsow scanner (GE Medical Systems, Milwaukee, Wisconsin). Details of the 
scanning parameters are shown in Table 1. 
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Table 1 . The MRI scanning procedure parameters. 

Scanning sequence Spin Echo (T1 weighted) 

Number of slices 54 

Pixel spacing 0.037x0.037 (cm-cm) 

Scanned region length (across the knee) 15.9 (cm) 

Slice thickness 2.5 (mm) 

Slice spacing 3 (mm) 



The Segmentation Procedure. A 3D tiled surface geometrical representation was 
generated using the software Amira (Indeed - Visual Concepts GmbH, Berlin, Ger- 
many), for the distal femur, the proximal tibia, and the insertion areas of the anterior 
(ACL) and posterior cruciate ligaments (PCL). 

A segmentation of the MRI data set was performed with an entirely manual 2D 
segmentation technique. For each slice, the outer contour of the structures of interest 
was detected and outlined, as shown in Fig. 1 . The resulting stacks of contours were 
interpolated to generate polygonal surfaces which represent the outer boundary of the 
objects to be modelled. The model used for the kinematic analysis is shown in Fig. 2. 

Kinematics. Series of lateral images were acquired at the frequency of 6 samples per 
second with a standard fluoroscope (SBS 1600, Philips Medical System Nederland 
B.V.). Images of a 3D cage of Plexiglas with 18 tantalum balls in known positions 
and of a rectangular grid of tin-leaded alloy balls 5 mm apart were collected in order 
to calculate respectively the position of the camera focus and the parameters neces- 
sary for image distortion correction. This was obtained using a global spatial warping 
technique[8]. An established technique for 3D kinematics analysis of a known object 
from a single view was implemented [9] (Fig. 3). Bone poses in space were obtained 
from each fluoroscopic image by an iterative procedure using a technique based on 
tangent condition between projection lines and model surface. Previous validation 
work on prosthesis components [9] showed that relative pose can be estimated with an 
accuracy better than 1.5 degrees and 1.5 mm. 




Fig. 1 . Outlined contours of femur, tibia, and ligaments in two slices of the MRI data set. 
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Fig. 2. Distal view of the femur (left) and proximal view of the tibia (right) model. The areas of 
insertion of ligaments are the regions on the femur and the tibia. 




Cruciate Ligament Geometrical Models. The geometrical models of the cruciate 
ligaments differ for dimension, shape of insertion and twist: 

1. Bi-dimensional - untwisted . The insertions were modelled as follows: a) the line 
that minimizes in a least squares sense the points of the insertion area was calcu- 
lated, b) insertion segment was identified on this line between the anterior and pos- 
terior limits of the insertion surface, and c) 25 uniformly distributed points were 
identified on this segment. Thus, 25 fibres were modelled for both ACL and PCL. 
In both ligaments the fibres connected these points of the insertions from the most 
posterior to the most anterior on both femur and tibia, i.e. the most posterior point 
of femur insertion was connected to the most posterior point of the tibia insertion. 

2. Bi-dimensional - twisted . The insertions and the fibre were modelled as in model 1 
except for the fact that the points from the most posterior to the most anterior of the 
femur were connected to those from the most anterior to the most posterior of the 
tibia. 
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3. Three-dimensional - rectangular insertions - untwisted . The insertions were mod- 
elled as follows: a) in the plane approximating the insertion points in a least square 
sense a rectangle including 80% of these points was estimated, and b) a 5x5 uni- 
form grid of points was identified on the rectangle. In both ACL and PCL the 25 
fibres connected points of the insertions with no twisting. 

4. Three-dimensional - rectangular insertions - twisted . The insertions and the fibre 
were modelled as in model 3 except for the fact that in both ligaments a twist angle 
of 90° was introduced. 

5. Three-dimensional - circular insertions - untwisted . The insertions were modelled 
as follows: a) in the plane approximating the insertion points in a least square sense 
a circle including 80% of these points was estimated, and b) a 25 uniformly dis- 
tributed points were identified on the circle. In both ACL and PCL the 25 fibres 
connected points of the insertions with no twisting. 

6. Three-dimensional - circular insertions - twisted . The insertions and the fibre were 
modelled as in model 5 except for the fact that in both ligaments a twist angle of 
90° was introduced. 

For each model, for each single fibre, the strain, £ , was calculated as follows: 



were L is the length of the fibre at time sample t , and L (} is the maximum length 
the fibre reached during the motor task. 



The modelled PCL always showed a larger elongation, with an average strain of about 
28% versus 16% of the ACL. The strain calculated for the fibre approximately con- 
necting the mean point of the insertions of the ACL and PCL was equivalent for all 
ligament models. The range of the strains calculated for the ACL [-14%;-18%] and 
for the PCL fibres [-20%;-37%] was similar for the two three-dimensional models, 
but was smaller, [-15. 9%;- 16.2%] for the ACL and [-27%;-30%] for the PCL, when 
calculated from the bi-dimensional model. The geometrical distribution of the strain 
over the ligament section resulted model-dependent. 

The strain calculated for the other fibres resulted also model dependent, in particu- 
lar the bi-dimensional models produced different results with respect to the three- 
dimensional ones. 

For the bi-dimensional models (Fig. 4) the PCL showed the largest strain at the an- 
terior fibres independently from the twist. The strain of the central fibres of the ACL 
was the smallest when untwisted and the largest when twisted. 

The strain behaviour of the fibres was similar for the two three-dimensional models 
(Fig. 5 and Fig. 6). The largest strain was observed for the ACL at the medial fibres 
when untwisted, and for the postero-lateral ones when twisted, for the PCL the medial 
fibres when untwisted and for the posterior fibres when twisted. 




( 1 ) 



Results 
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a b 

Fig. 4. The maximum value of strain over the section of modelled ligament during the execu- 
tion of the motor task for each of the 25 fibres is plotted for model 1(a) and model 2(b). 



ACL PCL ACL PCL 




a b 

Fig. 5. The maximum value of strain over the section of modelled ligament during the execu- 
tion of the motor task for each of the 25 fibres is plotted for model 3 (a) and model 4 (b). 



ACL PCL ACL PCL 




a b 

Fig. 6. The maximum value of strain over the section of modelled ligament during the execu- 
tion of the motor task for each of the 25 fibres is plotted for model 5 (a) and model 6 (b). 
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Discussion 

Six different cruciate ligament models were compared using parameters from a single 
selected subject analysed as accurately as possible. Plane, rectangular and circular 
sections were considered, and the mechanical effect of the anatomical twisting of the 
ligament fibres was also investigated. 

The strain range of the modelled fibres was dependent on the bi-dimensionality of 
three-dimensionality of the model, but was not relevantly influenced by the shape of 
three-dimensional model adopted. The more conventional bi-dimensional model [10] 
showed the largest differences from the two three-dimensional ones. No significant 
difference could be highlighted between the rectangular insertion and the circular 
insertion three-dimensional models. The twist showed significant influence in the 
strain distribution for each model. The described system, used to validate the cruciate 
ligament model, allows to reconstruct the in-vivo 3D kinematics of knee joint, com- 
bining 2D bio-images and subject-specific geometric parameters of the chosen bio- 
logical structures. The accuracy of this 3D reconstruction was evaluated previously 
for total knee replacement subjects. The best way to validate the method would be 
highly invasive for the subject under analysis. A phantom study would not add any- 
thing to the accuracy tests already performed. 

In conclusion, when only the mean magnitude of the fibres elongation is to be cal- 
culated the selected model does not considerably affect the results. Instead, the model 
should be accurately selected when the geometrical distribution of the strain over the 
section of the ligament is required, i.e. when the strain is used for the calculation of 
the load applied to the joint by the ligament [11]. In this case, a three-dimensional 
model is suggested, independently from the selected insertion shape, and the anatomi- 
cal twist of the fibres has to be taken into account, as it strongly influences the strain 
distribution over the section. This study on cruciate ligaments knee models is a start- 
ing point for the development of a useful teaching and clinical tool for the investiga- 
tion of the knee function in physiological conditions. 
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Abstract. Simulation of soft tissue behavior for surgical training sys- 
tems is a particularly demanding application of deformable modeling. 
Explicit integration methods on single mesh require small time step to 
maintain stability, but this produces slow convergence spatially through 
the object. In this paper, we propose a multigrid integration scheme to 
improve the stability and convergence of explicit integration. Our multi- 
grid method uses multiple unstructured independent meshes on the same 
object. It is shown that, with the proposed multigrid integration, both 
stability and convergence can be improved significantly over single level 
explicit integration. 



1 Introduction 

Simulation of the behavior of soft biological tissue for surgical training or plan- 
ning systems is a particularly demanding application of deformable modeling. 
Tissue is highly nonlinear, commonly showing an exponential stress-strain rela- 
tionship [1], Large deformations of 100% or more are possible [1]. Yet, simulation 
requires realistic behavior in real time. High accuracy at interactive speeds is 
necessary for intra-operative surgical planning (e.g., [2]). 

Recent research in soft tissue modeling has focused on finite element methods 
(FEM) for accuracy. Initial efforts (e.g., [3,4,5]) used linear models because of 
their speed and because they permit significant offline pre-computation to reduce 
runtime computation. Linear methods cannot handle large deformation or mate- 
rial nonlinearity, however. Cotin et al. approximated nonlinear elasticity within 
an otherwise linear model [6]. Picinbono et al. [7] and Muller et al. [8] extended 
linear methods to account for rotation that occurs with large deformation. 

In nonlinear FEM, strain is nonlinear in displacement (geometric nonlin- 
earity) while stress is nonlinear in strain (material nonlinearity). Szekely et al. 
[9] used a large scale multi-processor computer to obtain real-time performance 
while modeling both types of nonlinearity. Zlruang and Canny [10] modeled geo- 
metric nonlinearity for 3-D deformable objects, but not nonlinear material prop- 
erties. Zlruang used mass lumping to produce a diagonal mass matrix for speed. 
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Wu et al. [11] also used mass lumping, but incorporated both types of nonlin- 
earity. In this paper, we use the same modeling engine as [11]. 

Each of these nonlinear implementations uses explicit time integration for 
real-time response. In the animation community, implicit integration methods 
have become popular, spurred in part by the work of Baraff and Witkin [12]. 
While implicit methods permit large time steps in animation, convergence is 
typically too slow for real-time simulation of 15 or more frames per second if 
nonlinear behavior is accommodated. However, two recent papers have used 
implicit integration with finite element models for simulation. Muller et al. [8] 
warp a precomputed linear stiffness matrix according to the local rotation of the 
material to account for large deformations. However, their method only permits 
approximation of nonlinear stress-strain relationships, and cannot accommodate 
a general strain energy function. Capell et al. [13] use a multiresolution subdivi- 
sion framework that permits real time deformation, but this does not have the 
accuracy or flexibility of, e.g., tetrahedral meshes. 

Explicit integration schemes use much less computation per integration step, 
but require a small step size to maintain stability. The critical step size is propor- 
tional to element size and inversely proportional to the square root of stiffness 
[14]. Although the softness of deformable objects makes explicit integration fea- 
sible, there is a tradeoff between spatial resolution, i.e. the fineness of the grid, 
and the material stiffness that can be achieved. 

The other major disadvantage of explicit integration of dynamic models is 
that information is propagated spatially at only one layer per time step. Although 
the user’s motion while interacting with a surgical simulation is likely to be slow, 
the “slinky” -like motion that results from slow propagation is visually unrealistic. 
The stress concentration that results from local perturbations can also cause 
instability before the stress is distributed by the integrator. 

In this paper, we take a multiresolution approach to integration. Although 
there has been significant effort in multiresolution approaches in deformable 
modeling (e.g., [11,15,16]), past work has emphasized adaptation for providing 
local detail. Instead, our integration scheme maintains simultaneous multiple res- 
olutions of a triangle or tetrahedral mesh of the same object. The scheme is an 
extension of traditional multigrid methods in which stress due to a local pertur- 
bation is rapidly distributed through an object by restriction onto coarser mesh 
representations. After integration on the coarser meshes, results are projected 
back to finer meshes for refinement. The result is significantly faster convergence 
and better stability than single level explicit integration. 



1.1 Multigrid Iterative Solver Background 

Multigrid (MG) methods are also called multilevel, multi-scale, aggregation, and 
defect correction methods [17]. In the beginning, they were invented as one kind 
of iterative scheme to solve for X in 



T(X) = b 



(1) 
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where T can be either a linear or nonlinear operator mapping the current state 
X to internal stress, particularly for partial differential equations (PDE). In 
contrast to single level iterative solvers, MG uses coarse grids to do divide-and- 
conquer in two related senses. First, it obtains an initial solution from the fine 
grid by using a coarse mesh as an approximation. This is done by transforming 
the problem on the fine grid to the coarse grid and then improving it. The coarse 
grid is in turn approximated by a coarser grid, and so on recursively. The second 
way MG uses divide-and-conquer is in the frequency domain. This requires us 
to think of the error as a sum of eigenvectors of different frequencies. Then, 
intuitively, the work on a particular grid will attenuate the high frequency error 
by averaging the solution at each vertex with its neighbors. There are three 
operators used in the conventional MG iterative solver [18]: 

— The smoothing operator G(-) takes a problem and its approximated solution 
and computes an improved XW using a one-level iterative solver such 
as Gauss-Seidel or Jacobi methods. This smoothing is performed on all but 
the coarsest mesh, where an exact solution of (1) is formed: 



(i) _ /G(X«,6«),^(m-l) 
\T- 1 (6( Tn - 1 )),t= (to — 1) 



(2) 



The superscript (i) indicates G(-) operates on grid level i, for m mesh levels 
numbered from 0 (finest) to m — 1 (coarsest). 

— The restriction operator R(- ) takes the residual of (1), T(X^) — 2>W, and 
maps it to frb+i) on CO arser level i + 1. 

b (i+1) = R(T(X W ) - 6 (i) ) (3) 

— The interpolation operator P(-) projects the correction to an approximate 
solution XfW on the next finer mesh. 



X ii) = - P(X (i+1) ) (4) 

The three operators manage a series of grids with different resolutions and form 
a “V” shape algorithm as shown in Figure 1. 



1.2 Overview 

In this paper, we extend traditional multigrid methods to apply them to explicit 
integration of dynamic models of the form 

MU + DU + R{U) = V (5) 

where, for a finite element model, M and D are mass and damping matrices, 
respectively, U is the nodal displacement vector, R is the internal stress oper- 
ator, and V is the traction force vector. The methods operate on unstructured 
independent triangle or tetrahedron meshes. We developed and evaluated four 
variations, each using a single level explicit integrator in place of the smoother 
G(-) of traditional multigrid methods [19]. Different restriction and projection 
operators were used, producing different stability and convergence properties. 
Here we present the single method which combined the best properties of the 
four. 
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Finest Level: 




XO = G[T(X0)] 



R: restrict residual 
P: interpolate correction 
G: smoother 



Coarsest Level: 



i = 4 



X4 = T_inv(b4) 



Fig. 1. A V-cycle of traditional multigrid iterative solver solving equation (1). 

2 Methodology 

In the current implementation, we use the FEM engine of Wu et al. [11]. De- 
tails can be found in that reference. The FEM model uses a total Lagrangian 
formulation to accommodate geometric nonlinearity and multiple strain energy 
density functions have been implemented to represent the hyperelastic behavior 
of soft tissue. Mass lumping and Rayleigh damping are used to create a diagonal 
mass matrix M and associated damping matrix D. Consequently, the method 
scales computationally as O(n), where n is the number of nodes in the 2D or 3D 
mesh. 

2.1 Mesh Correspondence 

Our multigrid method uses unstructured triangle or tetrahedron meshes. The 
meshes are independent, i.e., no vertex needs to be shared between levels. The 
in levels of mesh, S A for * = 0, m— 1, are discretized versions of the continuous 
domain S. They can be chosen so that the numbers of vertices and elements 
in the levels form a geometric series with a predefined ratio. Because the FEM 
engine has linear computational cost O(n), the geometric series is the key to 
achieving 0(n o) cost as will be shown in section 2.3, where n o is the number of 
vertices in the finest mesh Off-the-shelf 2D and 3D mesh generators, such as 
Triangle [20] or Gmsh [21], can be used to triangulate or tessellate deformable 
objects, respectively. Both programs can control the output mesh density by 
setting a resolution parameter. 

The vertices and elements in each are mapped to other levels through 
an initial mapping process. Each vertex on level i has host elements on levels 
i — 1 and i + 1, for * £ [1 ,m — 2]. Each node in and S^ m_1) only has 
host elements in and g*™ -2 ), respectively. The correspondence is invariant 
throughout the simulation. As shown in the left of Figure 2, the host element of 
vertex P Q is identified as face {pi,P 2 ,Ps), which encloses Pq. The same is true 
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|P - cl| = |P - c2|, ml > m2 



ml/|P-cl| > m2/|P-c2| 



Fig. 2 . Left: Mesh correspondence in 2 D. Pi and p, are vertices in coarse mesh 
and fine mesh S^, respectively. The thick solid line and the thick dashed line in- 
dicate the boundaries of domains S*T and S (l_1 \ respectively. The fine face (pi,p 2 ,P3) 
is identified as the host element of coarse vertex Po. Fine vertex po is associated with 
coarse face (Pi,P 2 ,p3). P4 is located outside of the contour of fine mesh. Its host 
element is identified as face (p4,ps,P6), because it is the closest triangle to P4 mea- 
sured by ( 7 ). Right: Mapping boundary vertices, ci and C2 are the centroids of 
faces (pi,P2,P3) and (p2,P4,Ps), which are boundary elements. P is outside of the dis- 
cretized domain. Although P has the same distance to ci and C2, i.e. ypp = 7_p,2, P’s 
host element is classified as face fpi,P2,P3) because ■^ L - > -2i2_, 



for face (Pi, P 2 , P 3 ), which hosts vertex po- The weights for each coarse node are 
stored in a vector u> of size npe x 1, where npe stands for number of nodes per 
element. All interior nodes of coarse mesh S*T lie in the convex hull formed by 
the fine mesh S^ 1_1 ), and vice versa. Thus the associated weights for each coarse 
node ujj are non-negative and sum to 1. ujj is computed as 



Triangle mesh 
-yi, Tetrahedron mesh 



where Aj is the area of the sub-triangle involving the test point and two points 
from the host triangle. Vj is the volume of the sub-tetrahedron involving the 
test point and three points from the host pyramid. Thus the location of node 
pQ 1 can be represented as a weighted sum of the positions of corner nodes of its 
hosting element 7 on or those of its hosting element 77 on S^ 1+1 \ 




E npe 
0 = 1 

E npe 

0 = 1 



(i-l) 
7 ,j 






(i-l) 

1,0 



b+i) 

V,0 



UJ 



(»+l) 

V,0 



, de 7 G S**- 1 ) 
, ele v G S^ 1+1 ^ 



( 6 ) 



where p is the position vector of a node. 

Only boundary nodes of the coarse mesh are allowed to have negative weights 
or weights that are larger than 1. An element elei is the host element of the 
boundary node 8 on another mesh if elei encloses 9. If 9 is outside all the elements 
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on another mesh, then 9 is mapped to the boundary element elei on another level 
which maximizes 



nrii TOfc 

= max 

7 8,i k 7 8,k 



( 7 ) 



7 e,k = | \Pe 



1 

npe 



npe 



1=1 



7 $ t k is the Euclidean distance between 9 and the centroid of element k. As shown 
in the right half of Figure 2, if a vertex P is equidistant to the centroids of faces 
(.Pi 5 P 2 , Ps) and (P 2 ,P 4 ,P 5 ), face (pi,P 2 ,Ps) is identified as the host element if this 
face has a larger area. 

Each coarse vertex associates with a hosting element in the finer mesh. It is 
important that the hosting elements for the vertices in S ( b cover S^ 1-1 ). This 
guarantees that the state of every fine vertex can be represented in S*- 1 ) so that 
detail will not be lost during the restriction process R(-). This constrains the 
resolution ratio r between the number of vertices in and Sh -1 ). In practice, 
we have found r = 0.5 to be an effective ratio. 

Franklin’s pnpoly is used to test whether a point is inside of a triangle. The 
test of a point being inside a tetrahedron is conducted by checking if all four 
sub-tetralredra have positive volumes provided that the tested tetrahedron is 
not degenerate. After the correspondence between vertices and elements on dif- 
ferent levels is constructed, the uij’s are used as the weights when restricting or 
interpolating quantities between mesh levels. 



2.2 Multigrid Time Integrator 

The proposed Multigrid time integrator is called Restrict-and-Interpolate-State- 
or- Variation (RalSoV). Similar to the MG iterative solver, RalSoV uses a re- 
striction operator R!_ 1 (-) in the down stroke of the V-cycle in Figure 1 and an 
interpolation operator Pj^O) in the up stroke to transfer the problem from one 
level to another. Also in the up stroke, it applies a solution operator G(-) on 
each level mesh S*A. 



Restriction Operator. At current time step n, Rj^G) assigns either the state 
of each node or the variation of nodal state ZLYb -1 ) 

AX^-V = xg~V - X%ZP ( 8 ) 



in S^ 1 P to that in S 1 ^ in a weighted sum fashion. 

ni / w( i-1 h _ , ,(* — 1) yd -1 ) 

R ; _l(A 7 X -y,j > 



= 



| V • 0^7 1} AX (i ~ x) 



X. 



( i ) 



r; _ 1 (A4 i - 1) ) = x^ + E^ 7 , J 






(*— i) 



AX, 



J 70 

(i- 1) 

70 ’ 



70 

otherwise 



> e 



(9) 
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where the subscript 0 and (7, j) denote vertex index 0 on level i and the j th 
corner node of hosting element index 7 on level i — 1, respectively. The weight 
is computed as in section 2.1. The switch is conducted by checking the 
condition 

\Y,^j 1)AX ^ 1) \ >e ( 10 ) 

j 

When (10) is true, the weighted sum of the nodal state variation in the finer 
mesh is large. This indicates the regional nodes in Sb -1 ) are subjected to large 
stress, since the state differential X « l ar S e - R( - ) will linearly project 

those states into the corresponding vertices in the coarse mesh Sb). If (10) is 
false, the algorithm regards these finer nodes as converging to their steady state. 
By mapping the variation, RalSoV suppresses the high frequency component of 
the residual R(U) — V of (5). This improves the convergence rate of the system. 
For the proof of the linear problem please refer to Demmel [18]. 



Interpolation Operator. Similarly, the interpolation operator P] 1 (-) 
projects either Xb) or Z\Xb) in S (l ) to that on level i — 1 using weights 



p !_i (xP) = Zjv% x % 



b) v(i) 



x, 



b-1) 



|Z\X ( 



b-1) I 



> e 



Xg l 1) + Pj_ 1 (Z\X^) = Xg 1J + J2j u'rjAX'filj , otherwise 






di-i) 



b) 






(n) 

6 is the nodal index in S (l 1 h r is the hosting element index in Sb). The switch 
condition is 

" >e (12) 



l^l 



Although the formulations of (10) and (12) are different, their goals are consis- 
tent. By measuring the state variation in the finer mesh, we either project the 
state to relax the finer mesh local distortion when the nodal tangent is large, 
i.e. (12) is true, or interpolate the state variation to suppress the local residual 
in Sb -1 ) for faster convergence when (12) is false. For the rest of this paper, we 
will simplify the notation of R|_ 1 (-) and P^G) as R(-) and P(-), respectively. 



Solution Operator. RalSoV applies a conventional explicit ODE solver, such 
as forward Euler or Runge-Kutta, as the solution operator G(-) in each Sb). 



X 



(i) . 

n+1 



G(xb ) ,pb)) ; 



X 



b) = 

n — 



u 

u 



ib) 



(13) 



where pb) i s the right hand side of (5) at level i. The state of each vertex is 
a 6 x 1 vector including its displacement and velocity. This is different from 
the conventional MG iterative solver, which uses a direct solver at the coarsest 
level and a Jacobi-like pre-conditioner as a smoother in the other levels. The 
solution operator must be stable at each level. This condition ensures that XX 
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is closer to the steady state solution of (5), X& , than X„\ Currently, A is kept 
the same on every level for simplicity, though different step sizes could be used 
for different resolutions to accommodate individual stability requirements. If a 
nonlinear viscosity model is considered, the internal stress will require nonlinear 
mapping from both the position U and the velocity U in (5), i.e. R(U,U). 

The pseudo-code of RalSoV algorithm in the Appendix shows the method’s 
schematics. 

2.3 Computational Cost of MG Integrator 

Within each V-cycle, there are m time integrations, where m is the number of 
levels. The nonlinear FEM engine we use has a linear computational cost of 300 
FLOPS/node, for either 2D or 3D [11]. I.e., C(n ) = 300n where n is the number 
of nodes. This is much greater than the cost of the restriction and interpolation 
processes in the multigrid algorithms, 2 x 42 FLOPS/node for tetrahedra and 

2 x 30 FLOPS/node for triangle meshes. Thus, if the geometric series of the 
mesh nodes has ratio 0.5, the total cost for each V-cycle is roughly the sum of 
the FEM update cost at all levels, i.e. 

C Ra iSov{n (0) ) = c '( nW ) + Cost(mapping) 

( 2[C(n^) + 84n^] « 768n, 3D mesh (14) 

- \ 2[CV°)) + 60n (0) ] « 720n, 2D mesh 

The two sides of the above equation become equal as m approaches infinity. 
The increased computational cost is compensated by the stabler performance 
of multigrid integration, because in practice the multigrid algorithm can be in- 
tegrated with a time step larger than twice that of the single level integrator. 
This is shown in the 2D examples in the next section and the 3D example in the 
accompanying video. 

3 Evaluation 

In this section, RalSoV is compared to single-level integration method (SLIM). 
The better stability and the faster convergence of RalSoV are demonstrated 
with two examples. In the examples, four levels of triangle meshes are used to 
represent the discretization of the square membrane shown in Figure 3. The 
meshes were generated using Triangle [20]. SLIM evaluates only the finest mesh, 
So- The membrane has size 10 cm x 10 cm and all four sides fixed in space. 
It is modeled with a Mooney-Rivlin material [22] with constants Ci = 1.0 and 
Ci = 1.0. The damping ratio is 1.0. Both SLIM and RalSoV use an explicit 
fourth-order Runge-Kutta integrator to solve the system. 

3.1 Perturbation Response 

In the first example, with all four sides fixed in space and free of gravity, a central 
node of the membrane p'a'l is lifted instantaneously by 3.0 cm, and the node 
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Fig. 3. Four levels of 2D triangle meshes on a square membrane: So with 186 nodes 
and 321 faces, Si with 103 nodes and 165 faces, S2 with 55 nodes and 82 faces, and 
S3 with 31 nodes and 40 faces. Notice the numbers of faces at different levels form a 
geometric series. 

stays at the height throughout the simulation (Figure 4). This test simulates 
the step response of the membrane with concentrated large stress surrounding 
the lifted node initially. The initial displacement of paij causes the neighboring 
elements to be highly stretched. If a single level mesh is used to simulate this 
scenario, explicit time integration is going to be unstable unless a small step size 
A < A cr it = yj K c [14] is applied, where K ma x = is the largest 

eigenvalue of pseudo stiffness matrix. Note that, because of the hyper-elastic 
material property, K max increases as an element stretches. The time evolution 
of |T — Tod of the two algorithms is plotted in the left half of Figure 5, where 
T is the nodal stress tensor of caused by element deformation [11]. The 
stress of fixed boundary nodes is not included in the calculation. The single level 
algorithm is applied with time step 5 ms; RalSoV uses 10 ms as its time step. 

The initial lift is passed onto the nodes Pa} in coarser meshes S^ 1 ) to Si 3 ) that 
are closest to paj ■ Because the concentrated stress around pnj results in large 
variation at those surrounding nodes, RalSoV modifies those states directly. This 
largely regulates the stress concentration around Paj after one V-cycle, as can be 
seen in the left half of Figure 5. Eliminating the stress concentration eliminates 
the need to reduce the time step in order to maintain stability in the event of 
a large perturbation. Note that in Figure 5 SLIM does not distribute the stress 
on coarser levels as RalSoV does, and takes far longer for stress to reach steady 
state. In fact, if the perturbation were increased from 3.0 to 3.5 cm, SLIM would 
actually be unstable unless the step size were reduced further. 

After the strain distribution is smoothed, condition (10) is not satisfied for 
most nodes in the finer meshes. This “forgetting factor” eliminates the explicit 
mapping from the coarse meshes on successively finer levels. This improves the 
convergence of the result. 

3.2 Draping under Gravity 

In a second example, the membrane drapes under gravity with the four sides 
fixed. A time step of 10 ms is used for SLIM and RalSoV. The stress distribution 
in this case is smooth and this scenario is used to test the convergence rates of 
the two approaches. Figure 6 shows the initial and final states of the draping 
process. 
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Fig. 4. First: One node is lifted rapidly within one time step. Neighboring elements 
undergo very large deformation. Second: The displacement is restricted onto the coarse 
mesh, where the distortion is distributed over a larger region by the larger elements. 
Third: As the multigrid algorithm proceeds over multiple cycles, stress redistributes 
through the coarse mesh. Fourth: Effect of coarse mesh acts to spread influence over 
the fine mesh, eventually reducing its effect as steady state is approached. 

In this experiment, the gravitational force applies evenly on every element of 
the mesh so that strain distribution is also smooth. Condition (10) is not satisfied 
on all nodes throughout the simulation. RalSoV converges faster than SLIM as 
shown in right half of Figure 5. It combines the stability benefits of distributing 
large strain faster and the convergence attributes of conventional MG iterative 
solver while maintaining a computational cost bounded by 2C(n ^). 

3.3 Element Inversion 

Section 3.1 showed one way in which the multigrid approach can avoid instability 
by reducing the buildup of local stress. There is another way in which effective 
instability can occur. If an element is inverted by a perturbation large enough to 
cause a node to cross over the opposing edge (in a triangle element) or side (in a 
tetrahedron) within one time step, this will cause the model to go unstable. The 
multigrid method reduces this possibility in two ways. First, the distribution of 
stress distributes the effects of displacement over multiple elements, reducing the 
possibility that any of them will be inverted. Second, by restricting the effect onto 
a coarser grid with larger elements, the size of displacement necessary to invert 
an element is increased. These properties make multigrid more robust to this 
form of instability. This is demonstrated in a 3-D example in the accompanying 
video. 
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Fast 3.0 vertical lilt ot central vertex on a 10x10 membrane 




I 0 Single grid 
| » RSV 





Fig. 5. Left: One central node of the membrane is subjected to an initial vertical lift of 
3.0 cm. The results of RalSoV and SLIM are shown. Upper left: The time evolution of 
max | T — Tool, which occurs at the lifted node paj ■ Lower left: The average of \T — T^l 
over all internal nodes. Right: With four sides fixed in the space, the membrane drapes 
under gravitational force. Upper right: time trajectory of max\T — Too|. Lower right: 
mean\T — Too| versus time. 



4 Conclusion 

The multigrid algorithm RalSoV is easy to implement, provides clear benefits 
in stability and convergence over single level explicit integration, and provides 
real time performance. Multiple mesh levels of triangles or tetrahedra can be 
generated easily offline by off-the-shelf software. As indicated in equation(14), 
one MG cycle requires less than twice the cost of a single level method, so MG 
is cheaper to execute if it can stably simulate a dynamic system with the largest 
permissible time step max(hMG) > 2max(fe one _; e „ e (). Because of the improved 
stability of MG, this is likely to be the case. 

In our implementation so far, we have used the same integrator and time 
step at each level of the MG schemes. This is not required. Different integrators 
and step sizes have varying filtering properties that may produce advantages at 
different mesh resolutions. In fact, if the coarsest mesh has few enough nodes 
to permit implicit integration in real time, an intriguing possibility is to use 
explicit integration on the finer meshes to smooth the result while obtaining the 
stability advantages of implicit integration at the coarsest level. We are currently 
implementing an implicit solver for the FEM engine. 

The proposed method is not suitable for the analysis of the transient response, 
in which the details of the evolution of the state differential X are significant. 
However, it is reasonable to assume that the user’s movements in interacting with 
a surgical simulator will be relatively slow, so that dynamics can be neglected. 

A significant disadvantage of the multigrid method is that the mesh levels 
must be generated offline in practice. This prohibits, for example, local mesh 
refinement when tissue is cut. However, multigrid methods could be implemented 
within an overall hierarchical method such as the dynamic progressive mesh 
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Fig. 6. With four sides fixed in the space, the membrane drapes under gravity. Left: 
Membrane before gravity is applied. Right: Steady state after deformation due to grav- 
ity. 

scheme of Wu et al. [11]. Refinement could occur within the region of the cut. 
While a single level integrator would be necessary within that region, multilevel 
integration would still be possible elsewhere to maintain net benefits. 
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The RalSoV algorithm and the simulation video can be found at 

http : //www . medicals im. org/People/Xunlei/ symposium2004.html. 




A Suture Model for Surgical Simulation 



Julien Lenoir 1 , Philippe Meseure 2 , Laurent Grisoni 1 , and Christoplre Chaillou 1 

1 ALCOVE, INRIA Futurs, IRCICA-LIFL, UMR CNRS 8022, University of Lille 1 
59655 Villeneuve d’Ascq Cedex, France 
2 SIC, FRE CNRS 2731, Bat SP2MI, Bid Marie et Pierre Curie, BP 30179 
86962 Futuroscope Cedex, France 

Julien.LenoirSlif 1 . f r, 

Videos can be found on: http://www.lifl.fr/-lenoir 



Abstract. In this paper, we propose a surgical thread model in order 
for surgeons to practice a suturing task. We first model the thread as a 
spline animated by continuous mechanics. The suture is simulated via so- 
called “sliding point” constraints, which allow the spline to move freely 
while constrained to pass through specific piercing points. The direction 
of the spline at these points can also be imposed. Moreover, to enhance 
realism, an adapted model of friction is proposed, which allows the thread 
to remain fixed at the piercing point or slides through it. Our model 
yields to good results showing realistic behavior, robust computation 
and interactive rates. 



1 Introduction 

Suturing is a fundamental surgical gesture that any practitioner has to acquire 
and improve. This technique is useful for example during an ablation of an organ 
or for stitching a wound up. Besides, thanks to the growing power of comput- 
ers, we are able to offer interactive surgery simulations of realistic yet complex 
models. Consequently, it seems convenient to use surgical simulators to practice 
suturing in difficult contexts such as endoscopic surgeries [3]. We therefore want 
to design a dynamic surgical thread model for interactive simulation with which 
a surgeon can train the suturing. This model must be computed at 30Hz at least, 
for good visual effects. 

This paper is organized as follows: First, we discuss some previous work in the 
domain. Then, we explain our basis model and the new constraints we propose 
for suturing. Finally, we present some results before concluding. 

2 Previous Work 

Surgical simulation needs both models of organs and interactions. Thus, the tools 
that the practitioners use and their effects must be modeled too. However, most 
researches have concentrated on organs simulation only, whereas the modeling of 
certain tools remains an issue. Among these, the simulation of surgical threads 
has been studied only recently. The simplest approach to model a thread is to 
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use a mass/spring chain [10,4]. This technique is generally used in commercial 
simulators 1 . To avoid using a high stiffness for the stretching, some models even 
rely on a chain of rigid links [2], Pai [9] proposes a static simulation of a curve 
based on the cosserat theory. Moving a Frenet frame along the thin solid, they 
obtain a specific energy term measuring stretching and twisting deformation. 

Some algorithms have been proposed to handle suturing [4, 2] and knot tying 
[10] with such models. These however heavily rely on their discrete nature. For 
suturing, these methods do not result from any physical equations. At each time 
step, a point of the thread is linked to a point on the organ (considered as 
the piercing point). The concerned thread point can change to simulate sliding, 
but this modification is an arbitrary choice. Moreover, if the thread must pass 
through several piercing points, many thread points must be fixed and it is not 
clear if the resolution of the model remains stable. What is more, the overall 
movement results in “step-by-step” sliding. To hide this effect (visually and 
lraptically), the distance separating two successive points on the curve must be 
small, which induces a fine discretization and a penalized computation time. 
Moreover, this discretization loses the continuous property of the curve, which 
could have been kept by other approaches [14]. 

In a way similar to organ models, we want to design a model of threads, 
where all physical parameters remain continuous. Since during a suturing task, 
any point of the curve can potentially be constrained in the pierced holes, it is 
also desirable that the constraints should apply everywhere along the curve. For 
that purpose, we propose new types of constraints which allows the thread to 
slide with frictions through a point in a specific direction. The knot tying is out 
of the scope of this article and we only focus on sewing. 

3 Physical Simulation of Thread 

In this section, we briefly present our model of thread. However, more details 
can be found in [6] . We model the thread as a spline with few control points: 

n 

P(M) = ^qi(t)M s )- W 

i=l 

where s is the parametric abscissa, t the time, q; the control point positions, 
n the number of control points, and bi the basis functions specific to the spline 
type. We provide the curve with physical properties such as a continuous mass 
distribution and animate it using the Lagrange equation of motion. In this model, 
the coordinates q f with a £ (cc, y 1 z} of the control points of the curve are the 
degrees of freedom. 

One of the main interest of such a model, is that all the physical properties 
are defined in a continuous way. These include external forces or constraints, 
which can therefore be applied on any point of the curve and not only on the 
control points (which somehow do not lie on the curve) . The generalized form of 

See Surgical Science web site: http://surgical-science.com. 
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the forces applied at a point P are expressed as Qf = F and are thus auto- 
matically distributed among the control points. External forces include gravity, 
viscosity, etc. 

To control the deformations of the thread, we have to consider an internal 
energy that aims at physically structuring the spline. We can choose a set of 
springs regularly dispatched along the curve. The springs link two points of the 
curve (not necessarily control points) and are handled via the above-mentioned 
generalized forces computation. The use of various springs (including rotational 
springs) enable the control of stretching, bending and twisting. However, for the 
stretching energy only, a continuous approach exists [8] . 

To perform the physical constraints needed by a suturing task, we need a 
method which satisfies at most the constraints with no discontinuity in the phys- 
ical simulation. Instead of using unstable projection methods [11], we rely on the 
Lagrange multipliers method for its robustness. We obtain a linear system: 



/ M 00 ~Ll\ 
0 M 0 -LZ 
00 M -Ll 
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where L = ( L x L y L z ) is the constraint matrix, A the acceleration of the degrees 
of freedom, B a vector that sums the different contributions of all forces, E a 
vector coding the intensity of the violation of the different constraints. A are 
the Lagrange multipliers and each of them links a constraint to the degrees of 
freedom. The three symmetric matrices M form the generalized mass matrix M g 
and are computed in the following way (see [6]): 



V(i, j) G {l..n} x {l..n}, My = m. 



bi(s)bj(s) ds. 



(3) 



with m the mass of the spline. For n control points and c constraints, the overall 
system consists in 3n + c equations. 



4 Constraints for Suture Simulation 

We need several constraints to simulate the suture correctly. The first constraint 
is a sliding point constraint which allows our spline to pass through a specific 
point of an organ surface and to slide through it. Moreover, the thread direc- 
tion on the contact point is controlled to be orthogonal to the organ surface or 
directed by the needle inserted. This offers a more realistic simulation in which 
the thread does not turn freely around the contact point. To enhance realism, a 
local friction model is introduced on the contact point to reproduce both kinetic 
friction phenomena which brake the slipping and static friction sticking. 

4.1 Sliding Point Constraint 

The sliding point constraint allows the spline to slide through a specific point. We 
can consider that this constraint is similar to a fixed point one with a dynamic 
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abscissa parameter. It derives into three expressions (one for each coordinate) 
which we note g and are written as: 

g(q,q,M(t)) = P(s(t),t) -P 0 . (4) 



The equation g = 0 imposes that some point of the spline must be at the 
position Po, that we suppose fixed for now. Since it is a dynamic system, s(t) 
changes over time in order to be the right point of the curve that minimizes the 
energy of the constrained system. 

To avoid the drift due to numerical integration, we use the equations of the 
constraint g to formulate a second order differential equation. This provides a 
solution with a critical damping, known as a specific Baumgarte technique [13, 
1]. This leads to the constraint equations: 



jf 4+ i g = °- 



( 5 ) 



where At is the time step of the simulation. 

As s becomes a dynamic parameter, it also becomes a new unknown of the 
system which require a new equation. If we consider a perfect constraint, the La- 
grange theory imposes that the virtual power of the strain due to this constraint 
must be equal to zero. This is written as: 



<9g 

A -# =0 ' 

Os 



( 6 ) 



This theoretical framework is explained in more details in [12]. 

We decide to consider a different equation which gives us a direct relation 
between s and A to accelerate the resolution process. We allow that the effective 
work of the force generated by the Lagrange multipliers is not null and we 
represent it as an error. This error is due to an incorrect value of s and is 
applied to correct the dynamics of this parameter. This approach gives us an 
equation slightly different from Eq. 6: 



e .S + A.2i = 0. 

OS 



where the factor e is close to zero. The system becomes: 
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We decide to solve the system by decomposing the acceleration in two parts, 
one for tendency and another for correction: A = A* + A c [5] . The acceleration 
of tendency represents the acceleration without any constraint and the other 
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acceleration is the correction due to the constraints. This leads us to this new 
equation system: 

' MAf = B 1 
MA\ = B ' 

MAI = BZ 
' M g A c = L t A 
es = L T S \ 

L(A t + A c ) + L s s = E 

We replace the terms A c and s in the sixth equation by their expression 
respectively in the fourth equation and the fifth equation. These replacements 
yield an equation for A: 

' MA x t = B 1 
MAf = B y 
MA~ = 

< A C = M~ 1 L T X 
A 

es = Ls A 

, LM~ l L T X + ^-X = E — LA t 

To simulate a complete suturing task, the spline must pass through several 
piercing points Pq which implies several new variables s l . In practice, the algo- 
rithm works well and our system tolerates many sliding point constraints. The 
system is solved in 0(c.n 2 + c 2 .n+ c g .c 2 ) where c g is the number of unknown Sj. 

With only such constraints, the spline can freely slide through but also turn 
around all the piercing points without considering the organ surface. It appears 
necessary to avoid an inversion of the insertion type (that is, penetration in 
or exit out of the organ). For that purpose, it is convenient to design a new 
constraint which would insure the right direction of the piercing. 



4.2 Sliding Direction Constraint 

We need a specific constraint to impose at a sliding point that the spline is 
orthogonal to the surface of the organ. We thus define a direction constraint 
linked to a sliding point constraint. We just need the wanted vector direction 
To and the sliding point P(s(f),t) on which the direction is imposed. We first 
determine a local frame (P(s(t), t); T 0 , u, v) by computing u and v (vectors 
supporting the local plane tangent to the organ at P). We create a constraint 
that forces the direction of the sliding point to be orthogonal to u. Repeating 
this process on v, we constrain the direction of the sliding point to be in one 
direction orthogonal to u and v, it is thus in the direction Tq: 



<9P 

ci(q, q, t, s(t), u) = —(s (t), t). u = 0. 


(8) 


<9P 

c 2 (q,q,t,s(t),v) = —(s(t),t).v = 0. 


(9) 
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We constrain the direction of the tangent T(s,t) = but we let the 

intensity of this tangent free. Therefore, its norm will be set correctly by the 
energy minimization of the Lagrange equations. 

We still have a simulation of a thread that can freely slides through a point 
of an organ and in a specific direction. However, we now need to control the 
sliding. Frictions appear to be the most physically correct approach. 



4.3 Friction on Sliding Point Constraint 

At the point P(s,t), we compute the velocity of the point V = and its 
local tangente T = For any vector x (force or velocity), we can express its 
tangential component x t = (x.T)T and its normal component x ra = x — x # . 

To determine which friction model (static or kinetic) should be applied, we 
just have a check at the tangential velocity at the sliding point V t . 



• Static Case: If the velocity V t is below a given threshold, it is considered 
null, and we suppose that the frictions are static. In other words, the point is 
not supposed to move. Instead of computing the friction force which cancels 
the tangential forces (which would require a post-processing), we choose to 
replace the sliding point constraint by a fixed point one. The resolution of 
the system (cf. Eq. 8) gives the value of the Lagrange multipliers, which are 
related to the intensity of the forces which enforce the constraint. To respect 
the Coulomb friction model, we must check that the friction forces are in 
the friction cone. We compute both friction (At) and constraint (A n ) forces, 
and check if 1 1 At 1 1 < Ms 1 1 Anil, where Ms is the static friction coefficient. If 
the test succeeds, the static friction hypothesis holds and the computation is 
over. Otherwise, the frictions are pseudo-static. In that case, the constraint 
is replaced by a sliding point constraint. We keep the normal component \ n 
of the constraint force, and add a tangential force to the system: 

A friction — Ms 1 1 1 1 T7T 77 

II All 






Dynamic Case: If the velocity V) is not null, the kinetic friction model is 
immediately applied. The system defined with the sliding point constraint 
without friction is solved. The Lagrange multipliers give the normal force 
X n . We then compute the kinetic friction force, based on the kinetic friction 
coefficient M/M 



Ffriction — Mfc 1 1 



Vt 

'iivtir 



This algorithm is applied to all the piercing points. We have been able to 
check if static, pseudo-static or kinetic frictions were to apply at each point. If 
frictions are static everywhere, the computation of the forces is over. If however 
one or more points are in the pseudo-static or kinetic case, the computed friction 
forces have to be injected into the system which is solved again (to allow the 
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friction forces to be dispatched to the degrees of freedom of the curve). Thus, 
our algorithm generally requires two solving passes of the system of Eq. 8. 

It can be noticed that, in the kinetic case, we start the computation step 
by considering a sliding point constraint which is biased by e. Theoretically, 
the Lagrange multipliers give us a force orthogonal to the local tangent (i.e. 
\ t = 0) due to Eq. 6. However, by considering Eq. 7, we allow the constraint to 
work. We thus compute the effective normal force by projection and the non-null 
tangential force is just considered as a supplementary force applied to the thread 
(like deformation forces, collisions...). 



4.4 Interaction between the Thread and the Organ 

Since the beginning of this section, we have always considered that the points 
Pq of the organ were static. To take the motion and deformation of the organ 
into account, a quasi-static approach is possible. We consider Pg as constant for 
a given time step (this assumption is reasonable, since the motions are small). 
After all the sliding points Pg have moved, the sliding point constraints may be 
violated, but the stabilization scheme of Eq. 5 attracts the curve toward all the 
points. For the suture to modify the movement and deformation of the organ, 
we use the action/reaction law. At each Pg, the normal and friction forces are 
inverted and applied to the organ. 

A more precise and robust approach consists in assembling the thread and 
the organ in a same system [5]. (Pq,u®,v*) are then considered as variables of 
the organ. The constraint equations Eq. 4 and Eq. 8 can then differentiate the 
expression of (Pg, u®, v®) according to the motion of the organ. 

5 Results 

Our implementation takes place in a framework of surgical simulators [7] offering 
tools for collision process and self collision detection. This platform also gives 
access to a large choice of numerical integration methods which are necessary to 
compute the motion of our model. 

The simulations were performed on a 1.7GHz Pentium IV processor with 
512 MB of memory and an implicit Euler numerical integration was used. The 
test (figure 1) simulates a suture controlled by a needle. We use a spline with 
20 control points and two piercing points. The first figure presents the initial 
situation and the second one shows the simulation at some time later when the 
user pulled the needle. We can see that the thread has slid through the two 
piercing points and it is directed normally to the organ surface. The simulation 
takes about 30 ms which implies a simulation frequency of 34Hz. 

The spline interpolation allows the use of a limited number of control points. 
However, during a suture, it is important to give enough degrees of freedom for 
the curve to deal with all the constraints. We are currently working on a multi- 
resolution method which can not be described here due to space limitation. This 
method can locally increase the number of control points while removing others 
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Fig. 1 . Suturing task with a thread and a needle. The figure on the left shows a thread 
sliding along two piercing points. The figure on the right shows the state of the thread 
after the user has pulled it. 

in low curvature areas. In our first experiment, we see that, even for complex 
geometric configurations (such as a knot), very few degrees of freedom are added 
and the number of control points remains quite the same. This shows that 20 
control points is in practice sufficient. 

6 Conclusion and Perspectives 

We propose in this paper a dynamic model of surgical threads with specific 
constraints permitting the suturing task on objects like organs. These constraints 
are enabled to force the thread to pass through a sliding point, to impose the 
direction of the thread at this point, and to integrate frictions between the thread 
and the organ. 

Nevertheless the complete suturing simulation still requires a lot of work. 
First, it is essential to identify our model to the properties of existing threads 
to get a realistic simulation. For such purpose, the adaptation of the Cosserat 
energy seems inescapable [9]. Besides, the knot tying was out of the scope of this 
article. It requires the multi-resolution process that we have mentioned above. 
It is still under development but show promising results. In this context, we are 
also working on the cutting of the thread. We are confident that our thread 
model will become a powerful tool for future simulators. 
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Abstract. Surgical simulations with the aid of computers is a topic 
of increasingly extensive research. Real-time response and interactivity 
are crucial components of any such system and a major effort has been 
invested in finding ways to improve the performance of existing systems. 
In this paper, we propose the use of a new variant of Free Form Deforma- 
tions that supports discontinuities (DFFD) in complex real time surgical 
simulations. The proposed scheme can benefit from an a priori or, al- 
ternatively, interactive low resolution, physical simulation of the incision 
procedure that is encoded into the DFFD. The DFFD is then applied 
to the geometry of the tissue at hand, a two-dimensional skin shape or 
a three-dimensional volumetric representation, capturing the incision’s 
shape as well as the behavior of the neighborhood geometry over time. 



1 Introduction 

Surgical simulators have become an active research subject in recent years. Sur- 
gical simulators allow physicians to improve their skills inside a virtual envi- 
ronment prior to entering a real operating room. Such pre-operative training 
procedures have been shown to significantly improve the results of actual proce- 
dures [1], 

In order to maximize the potential gain in such pre-operative training, a sur- 
gical simulator should replicate the surgical environment as closely as possible. 
Conveying a realistic impression of even the simplest procedures is a difficult 
task. Because of its complexity, it should be broken into much smaller undertak- 
ings. One of the most important parts of any surgical simulator is to realistically 
and interactively animate the way tissue behaves under cutting operations. A 
virtual cutting module should supply, as a minimum, the following basic ca- 
pabilities. First, it should implement some mechanism for collision detection. 
Such a mechanism should control the location, direction and orientation of a 
virtual scalpel. Second, a cutting module should implement geometrical oper- 
ations that would progressively cut through the virtual model, modifying the 
topology around the cut and constructing new geometry as needed. Third, the 
cut model should present as accurately as possible physical behavior, including 
also capturing tissue behavior over time. 

The task of tissue cutting can be divided into two sub-tasks. First, the sur- 
face of the model should be split along the route of the virtual scalpel. Second, 
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the geometry in the vicinity of the cut should change, reflecting the shape of 
the cutting tool and the tension in the tissue. In this work, we suggest a frame- 
work that would combine these two subtasks into a single, unified operation. 
The method is based upon an augmented variant of Free Form Deformation [2] 
(FFD), which can admit discontinuities and openings into geometric models. 
A discontinuous FFD (DFFD) is continuous everywhere except at the incision, 
and hence it has the potential advantage of being able to continuously deform 
the geometry around the cut. Moreover, powerful analysis models such as Finite 
Element Methods (FEM) can be used to approximate the tissue’s response to 
the cutting operation and this response can be approximated and coded once 
into the DFFD. Alternatively, a physical simulation could be applied to a low 
resolution model, encoded into the DFFD during the interaction, and immedi- 
ately applied to the fine resolution representation of the geometry. Either way, 
the encoded DFFD captures the time-response of the tissue around the cut and 
etches this time-response to the neighborhood of the cut, as time progresses. 

While it can benefit from physical analysis such as FEM, the proposed 
method is purely geometric with low computational complexity. Consequently, 
the algorithm is capable of handling complex geometric models in real-time. 
However, since it, potentially, performs no direct physical simulation as part of 
the cutting operation, care must be taken in the construction of the deformation 
function in order to achieve convincing results. 

The rest of this work is organized as follows. In Section 2 we give an overview 
of the previous work on the problem of model cutting. In Section 3 we describe 
the proposed cutting approach; Section 4 presents a few examples. Finally, we 
conclude in Section 5. 

2 Related Work 

Throughout the years, the problem of cutting through geometric models has been 
tackled from multiple directions. In this section, and due to space constraints, we 
only consider a small subset of the relevant work. Some efforts considered surface- 
based meshes to represent the cut object. Ganovelli and O’Sullivan [3] suggested 
cutting surface-based meshes while re-meshing the cut to achieve the required 
level of smoothness in the cut. Since such splitting operations could degrade the 
quality of the mesh, [3] suggested the use of edge-collapsing to remove low-quality 
triangles from the mesh. A different approach aiming to preserve mesh quality 
under cutting operations was proposed by Neinhuys and Van der Stappen [4] . 

[4] suggested using a local Delauny-based triangulation step as part of the 
mesh-cutting process. Edge-flip operations are used on the faces affected by 
the cutting operation. The result of employing edge-flips is the elimination of 
triangles with large circumferences. The above methods operate on a polygonal 
approximation of the true free- form surface. A different approach, by Ellens and 
Cohen [5], directly incorporates arbitrary-shaped cuts into tensor product 13- 
spline surfaces. This approach has the advantage of operating over an inherently 
smooth surface, but requires some modifications to fit the standard definition of 
trimmed B-spline surfaces. 
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Other works considered a volumetric data model, usually in the form of a 
tetrahedral mesh for model representation. Using a volumetric data model is 
beneficial since it can represent both the outer surface of the model as well as 
its inner parts. Bielser et al. [6] described cutting through tetrahedral meshes 
based on the observation that there are only five different ways to cut a tetra- 
hedron topologically. A collision detection algorithm is used to find the position 
and direction of a cutting scalpel and map each cut to one of five available cut- 
ting configurations, adding additional tetralredra into the model and degrading 
the performance of the system over time. To prevent this problem, Neinlruys et 
al. [7] proposed locally aligning the mesh to the route of the virtual scalpel. The 
movements of the scalpel inside a triangle are recorded and the vertices adja- 
cent to the motion-curve snap onto it. Then, the triangles are separated along 
the aligned edges. Another effort, aiming at minimizing the amount of added 
primitives in a tetrahedral mesh during cuts, by Mor and Kanade [8], explored 
different subdivision strategies of a tetrahedral mesh. Model deformation was 
achieved by employing linear FEM over the cut model. Another problem that 
was tackled in [8] is the introduction of progressive cutting to prevents delays 
during the cutting procedure. 

A different approach, proposed by Forest et al. [9], treats cutting through 
tetrahedral meshes as a material removal problem. In [9] , tetralredra are removed 
from the mesh when hit by a pointing device. The method trivially conserves 
the tlrree-manifoldness of the tetrahedral mesh but also results in a loss of mass 
to the volumetric model. Since the fineness of the cut is tightly coupled to the 
fineness of the model, the result could present sharp edges that are less commonly 
found when cutting human tissue. 



3 The Algorithm 

The proposed algorithm operates in two phases. First, a discontinuous deforma- 
tion function is constructed and positioned over a local region of the surface of 
the model that is undergoing an incision procedure. In Section 3.1, we describe 
this process of constructing the deformation function. Following this, the geom- 
etry that is embedded in the spatial domain of the deformation function must 
be split and the deformation process over the cut geometry is performed over 
time. Section 3.2 describes the process of splitting the mesh as a response to the 
movement of the virtual scalpel; while Section 3.3 elaborates on the deformation 
process over time. 

3.1 The Deformation Model 

We briefly describe the deformation model that is used. Consider a trivariate 
function, F(u,v,w) : 1R 3 — » IR 3 , which maps a box-shaped parametric domain 
into a contorted box in Euclidean space. Deformations to any geometric model 
could be specified by properly designing the shape of the deformation function 
F. This deformation paradigm is known as FFD [2] and has been successfully 
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applied in many domains such as Computer-Aided Geometric Design (CAGD) 
and registration applications in medical imaging [10-12], to name a few. In this 
work, we used a tensor product trivariate B-spline function for the deformation 
function: 



l m n 

F(u, «,«,) = £££ PijkB°{u)B°(v)B°(w), (1) 

•£=0 j — 0 k — 0 

(u, V, w) G [Ujnin , Umax) X [VrniniVmax') X \Wmiri')W^max) 

where Pijk are the control points of F, and B{u) are the univariate B-spline 
basis functions, of order o, in all three directions. 

In the context of a cutting application, regular FFD suffers from a major 
shortcoming due to its continuous deformation nature. Being continuous, inci- 
sions and cuts cannot be introduced by the regular FFD paradigm. In [13], we 
circumvent this limitation by allowing the introduction of discontinuities into 
the FFD. We denoted this variant of FFD as DFFD for Discontinuous FFD. A 
DFFD function F(u,v,w ) will possess C -1 discontinuities along prescribed iso- 
parametric direction(s). To that end, a knot insertion procedure is used to insert 
d knots into the parametric domain of F yielding a potential C _1 discontinuity 
along the route of the virtual scalpel. As a result, two planes of adjacent con- 
trol polygons are added into the control lattice of F. Manipulating these points 
would let us model the shape of the discontinuity. The DFFD can be continu- 
ous at times and discontinuous at others, as needed. Hence, the insertion of an 
incision could be translated into a discontinuity in the shape that is formed by 
a cutting operation; see Figure 1. 




Fig. 1 . The application of a DFFD (b) to a polygonal mesh (a) could yield the incision 
insertion shown in (c). 

A DFFD can be applied to local parts of a geometric model, deforming it 
while incorporating discontinuities and openings into the geometry of the model. 
As the user traces a route on the screen with a pointing device or a virtual 
scalpel, a collision detection procedure projects each point and reconstructs its 
three-dimensional hit position on the surface of the geometry. The available 
hit points are used to construct a B-spline curve, C (u) , which approximates the 
movement of the virtual scalpel over the surface of the geometric model. In other 
words, curve C(v) prescribes the motion of the virtual scalpel over the surface. 
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Let T(v) = dC Ju ' > be the tangent vector of C(v) and let 0(v) be the orientation 
of the scalpel. Then, N(v) = 0(v) — prescribes the orientation of the 

scalpel in the normal plane of C(v). N(v) is typically normal, or close to normal, 
to the surface that undergoes cutting. Figure 2 presents these curves. 



N{v) 





Fig. 2. The DFFD is constructed with the aid of the scalpel’s path, C(v), and the 
scalpel’s orientation, 0(v), projected onto the normal plane of C(v). 



C(v), T(v) and N(v) fully prescribe the immediate behavior of the cut inser- 
tion. From these three vector fields the DFFD is reconstructed to follow the path 
of the incision. Nevertheless, the DFFD can also provide an emulation of the be- 
havior of the neighborhood of the cut over time due to internal tissue forces. 
Based on the type of tissue in hand, this time-response could be combined with 
these three vector fields, in order to encode a complete DFFD function. This 
DFFD not only captures the incision’s proper position and orientation but also 
the behavior of the neighborhood of the cut over time. 

The initial construction of the DFFD function F is based on C(v), T(v) 
and N(v). The control volume of F is built around samples along C(v) and 
the uniform distribution of the Pijk control points at each sample C t along the 
plane spanned by Nj , the sample on N(y), and T,j x iV,; , and is explained in 
[13]. Subsequently, the time response can be encoded by applying the a priori 
computed time response behavior coercing the cut to open and follow the pre- 
computed behavior of the tissue. 

Denote by T r the amount of time it takes for a cut tissue under tension to 
reach a steady state. The anticipated behavior of tissue could be encoded with 
a DFFD that starts at the position of the scalpel and goes back through all 
the events that the virtual scalpel captured after time T c — T r , where T c is the 
current time. Hence, the DFFD function F is of finite support, regardless of 
the complexity of the underlying geometry. Moreover, such a treatment of past 
events also captures the anticipated behavior of slow vs. fast moving operations. 
When the motion is slow, the cut will be fully open close to the scalpel, having a 
short window along C(v) that is governed by F. Alternatively, when the scalpel 
is fast moving, the cut will complete its opening far behind the current scalpel 
position, having all the events from time T c — T r all the way to time T c in a long 
window along C(v) that is governed by F. 
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3.2 Splitting the Geometry 

While the DFFD could be applied to surface or volumetric data sets, polygonal 
or polynomial, in this work we concentrate on triangular meshes. When cutting 
through a triangular mesh, the basic operation one needs to consider is triangle 
splitting. Any triangle that has vertices on both sides of the cut must be split 
before the DFFD mapping is applied to it. Moreover, vertices on the cut must 
be treated with care, moving them to the proper side of the cut. This splitting 
operation is conducted incrementally as the virtual scalpel moves, one triangle 
at a time and following the path of the curve along the skin surface. 




Fig. 3. The process of splitting the triangle and adding the deepness of the cut’s 
geometry is presented. The triangle in (a) is split into three triangles along the scalpel 
path (6), interior triangles are added in (c) and the DFFD is applied in ( d ). 



Splitting the triangle creates three new triangles; one triangle on one side of the 
cut, and two on the other side. These three new triangles replace the old triangle, 
which is then purged. Every triangle has both entry and exit cut locations. 
In addition, we examine how close the linear approximation at C(v) is. If the 
straight line between the entry and exit locations is sufficiently close to C( v), 
that triangle is split. Otherwise, the triangle must be refined recursively into 
smaller triangles before the split can take place. Figures 3(a) to (d) show an 
example of one such split with the new faces in red. 

In this application, we only process a skin surface. Nevertheless, we also seek 
to model the deepness of the cut, and model this new geometry on the fly. Toward 
this end, the entry and exit locations are duplicated and moved along the skin, 
in a direction of —N(v) and into the body. Triangles are then used to tessellate 
these interior vertices, all the way to the bottom of the cut. This process can be 
seen in Figures 3(c) and (d). 

Every vertex in the deepness of the cut is associated with texture coordinates, 
according to the distance from the beginning of the cut, along C(v) and its depth, 
along N(v). The typical texture that corresponds to this tissue’s cross section is 
then displayed on the side walls of the cut. 



3.3 Deforming the Geometry 

The DFFD is being updated incrementally as time progresses, clipping tail events 
that are below time T c — T r and adding new events that come in from the virtual 
scalpel. At each time step, all vertices of the model that are inside the current 
window of F are evaluated and deformed to their new position. As a result, the 
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total motion that each vertex undergoes is the consequence of all these small 
motions, in each time step, integrated over time as the window of the DFFD 
advances. All the vertices that are on the right side of the scalpel, including the 
vertices that are on the cut, which are part of the triangles on the right side, are 
mapped to the right, and likewise for the left side. 

In order to apply the DFFD, we need to parameterize the vertices of the 
skin surface in the window of F. In other words, we need to assign ( u,v,w ) 
coordinates to each and every vertex we intend to move. For every skin vertex 
Vi, let v = argmin ||C(t;) — Vi\\. Then, u and w are set by their distance to C(v) 

V 

along N(v) and T(v) x N(v), respectively. The ( u,v,w ) values are computed 
once for each vertex Vi, which requires proper parametrization of the DFFD 
function in the moving, V , direction. Since only the parametrization of new 
vertices, the ones that enter the moving window, need to be computed for every 
DFFD execution, the procedure is even less time consuming. Vertex, Vi, on the 
cut will move an amount that is the integral of all the small motions assigned 
to it via the evaluation of F, while V, is in the window of F. 

Now, using the DFFD, we calculate the new position of the vertex, and 
update the mesh. Vertices that lie outside the window, i.e. have u, v or w values 
outside of the domain, are not treated. This means that for proper continuity 
away from the cut, the mapping of F should preserve the vertices’ position and 
orientation along the boundaries of F. 

4 Results 

Figures 4(a) to (/) show six cutting profiles that may describe different tissue 
responses to the use of varying cutting tools. Figures 4(a) and ( b ) show cuts 
that leave wider gaps in the upper section of the incision. Figure 4(c) models 
a minimal cut at the level of the skin with a deeper cut in the inner tissue. 
Figures 4(d), (e) and (/) show different symmetrical and asymmetrical cuts. 
These cutting profiles demonstrate the flexibility with which DFFD can model 
cutting tools of varying geometry. 




Fig. 4. The DFFD framework is used to generate cuts of varying shapes with ease. 
The different shapes of the cut are encoded into the shape of the DFFD function. 
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Fig. 5. Several incisions that were produced using the DFFD algorithm, (a), (b) and 
(c) show cuts that intend to simulate plastic surgery procedures. 



Figure 5(a) shows a simulation of a cut that is used for facial plastic surgery, 
rhytidectomy. The skin around the ear is cut in this facial plastic surgery proce- 
dure. Figure 5(b) shows a different type of cut in the nostrils. Figure 5(c) mimics 
a cut from an eyebrow surgery. The examples were constructed to resemble im- 
ages of real facial surgerical procedures as shown in [14] . 




(a) (b) (c) 



Fig. 6. Models of human internal organs, (a) shows a stomach. ( b ) a liver and (c) shows 
a heart. An incision is introduced into the surface of all the models. 

In Figure 6(a), a model of a human stomach is shown. The cutting algorithm 
was used to simulate an incision in the virtual model. Figure 6(b) shows a model 
of a human liver, while Figure 6(c) shows a model of a human heart. In all 
these models, the DFFD cutting algorithm was used to incorporate cuts into 
the model. All three models are colored in shades of red. The coloring scheme is 
artificiality generated and does not intend to reconstruct the true colors of an 
original human organ. 

5 Conclusions and Future Work 

We have proposed an algorithm for real-time incision simulation that is suitable 
for both meshed surfaces and volumetric models, polygonal or spline-based. To 
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that end, free form deformations that support discontinuities, or DFFD, are 
employed. The use of DFFD lets us capture both the immediate geometry change 
due to the incision and the time-response of the tissue in the vicinity of the cut 
due to tension and internal forces. 

Although we implemented only the meshed surface option, the versatility of 
the DFFD allows the handling of volumetric meshes as well. The technique is 
simple and fast since no physical considerations are directly involved, allowing for 
real-time interaction. Future work could include the extension of the algorithm to 
support volumetric tetrahedral meshes and/or volumetric data sets. Volumetric 
representation could also improve the texture mapping capabilities over the walls 
of the cut, using the now readily available volumetric texture information. 

In addition, the incorporation of physical analysis into the DFFD during the 
interaction and the use of low resolution FEM analysis should also be explored. 
An incision into a flexible tissue that possess internal tension would be simulated 
by a DFFD with a large opening, rigid tissues could be represented by a small 
opening in the deformation function. One may record the time response of a 
tissue to an incision operation. Then, the time response could be fitted by a 
DFFD function, resulting in an approximation to the tissue’s response. Such a 
modelling task can be solved by any data-fitting scheme. Alternatively, one could 
use a low resolution FE three-dimensional grid to compute an approximated 
solution, then use the nodes of the grid as the control point for the DFFD 
function that can be applied to the high-resolution mesh of the model. In places 
where an exact tissue response is less important, such a combination could supply 
a low-cost and interactive solution. 
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Abstract. Interactively simulating nonlinear deformable human organs 
for surgical training and planning purposes demands high computational 
power which lacks in single processor machine. We build an interactive 
deformable objects simulator on a highly scalable computer cluster using 
nonlinear FEM and the novel multigrid explicit ODE solver which is 
stabler than single level schemes. The system consists of a graphical front 
end client on a workstation connected to a parallel simulation server that 
runs on a Linux cluster. After discussing the methodology in detail, the 
analysis of the speedup and preliminary results are presented. 



1 Introduction 

Advances in networking, visualization, and parallel computing have started to 
replace batch mode processing for computationally intensive applications [3] and 
to allow the real-time visualization of large data set [1,2,4]. Soft human organs 
are often subject to large nonlinear deformation during vast surgical procedures. 
Accurate and interactive simulation of such complex behavior in high resolution 
is a new challenging application which requires a large volume of computation. 

To provide the necessary computational power for the real-time solution of 
large scale dynamic systems, Szekely et. al. [6] have constructed a cluster of 
machines to perform endoscopic surgical simulation based on FEM. The cluster 
with predefined 3D lattice topology requires the object geometry is manually 
partitioned in a similar manner. The inter-processor communication is also pre- 
designed to 6 neighbors. 

1.1 Contribution 

To overcome the limitation and make the high performance surgical simulator 
design more flexible, we present a scalable deformable object simulation proto- 
type incorporating a real-time nonlinear FEM engine and multigrid (MG) time 
integration scheme. 



S. Cotin and D. Metaxas (Eds.): ISMS 2004, LNCS 3078, pp. 124—133, 2004. 
(c) Springer- Verlag Berlin Heidelberg 2004 
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Highly Scalable and Portable Implementation. This system is imple- 
mented in C++ on an existing cluster system, Millennium,, developed and con- 
structed at University of California, Berkeley 1 . In Millennium, each comput- 
ing unit, or so called node, is a Pentiumlll Linux workstation holding its own 
CPU(s) and memory. High speed network, e.g. MyriNet 2 , links each node to a 
network switch. Topologically, it has virtually a star shape where every node 
connects with every other node. This dramatically enhances the flexibility in 
data distribution and load balancing over 3D lattice topology designed in [6]. 
Inter-processor communication is managed by Message Passing Interface (MPI) 3 
which is portable through many parallel platforms and is binded to Fortran77 
and ANSI C for broad compatibility. 



Stabler Explicit ODE Solver. In section 4.1, a novel multigrid explicit ODE 
solver utilizes multiple meshes with different resolutions to render the object and 
simulates the deformation. The method achieves faster information propagation 
than single level explicit ODE solver and thus improves simulation stability. 



Remote Visualization and Interaction. With minimized bandwidth con- 
sumption, the data transmission between the visualization client and the com- 
putational server can be handled by widely available 100Mbps Ethernet. This 
will enable further research into remote scientific education and interactive large 
data set visualization addition to tele-surgical training/planning applications. 

2 System Design 

Our system consists of a front-end client on a graphical workstation connected 
to a cluster simulation server that runs on a Linux cluster as shown in Fig. 1. 

2.1 Simulation Server 

The simulation server implements the data manager and the actual simulation 
engine on a master node Pq and a set of slave nodes Pi, ...Pjv- The data manager 
on P 0 reads in the deformable object geometry, the constraints, and material 
parameters. It then launches a server thread to communicate with the front-end 
clients and to serve requests from slave nodes and front-end client during each 
time step. Each slave node Pi, (i > 0) runs in a compute- and- communicate loop. 
During the compute phase P, performs either single level ODE solver when only 
one level is simulated or MG time integration when multiple meshes are engaged. 
During the communication phase Pi first sends/receives the ghost nodes detailed 
in Sect. 3.1 to/from the neighboring nodes, and then Pi sends its hosted vertices’ 
new state to Pq. 

1 http : //www. millennium. berkeley . edu 

2 http : // www . myr i . com 

3 http : //www-unix .mcs . anl . gov/mpi/ 
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Millennium 




Socket 

TCP/IP 




OpenGL UI 



Fig. 1. The system diagram of our interactive parallel deformable object simulator. The 
workload is distributed in Pi ... Pn and the result is gathered by a central processor 
Pq. The inter-processor communication is conducted by MPI. 

2.2 Graphical Front-End Client 

The graphical front-end is the user interface to the simulator for both visualiza- 
tion of organ deformation and interpretating user interaction with the organ by 
mouse. The client receives the data stream from Pq. The received data contains 
only the vertex positions of size N x DOF x sizeof (float), because the vertex 
connectivity information is initially transmitted to minimize network congestion. 
For example, simulating the largest mesh with 12,978 elements and 6,657 ver- 
tices at 25fps requires 2MByte/sec bandwidth, which can be handled by 100Mbps 
Ethernet easily. This means that the visualization client can be remotely located 
away from the computation cluster which broadens the application. 

3 Parallelization 

Evaluating element’s internal stress T [9] is much more costly than explicit in- 
tegration on each node. METIS 4 is used to convert mesh’s element connectivity 
into a graph and then evenly partition the graph by assigning each element 
unique rank indicating which processor it belongs to. The mesh partitioning is 
done in every mesh level so that each slave CPU posesses one slice of the mesh 
at each level. 

3.1 Load Balancing and Distributing Ghost Nodes 

Information from neighboring portions of the grid(mesh) “owned” by other pro- 
cessors must be communicated to the given processor; this communication is 
usually done through the concept of ghost nodes [5] as shown in Fig. 2. E m 
and E n are two elements with rank p and q , respectively. Since V 3 has rank q 

http : //www-users . cs .umn . edu/'karypis/metis/ 



4 





An Interactive Parallel Multigrid FEM Simulator 



127 




Fig. 2. Left half: Identify ghost nodes. Right half: Distribute ghost nodes. 

which differs from that of E m , V 3 is a ghost node on processor q. Similarly, V\ 
is a ghost node on CPU p. If the nodal rank is simply assigned according to its 
hosting element’s rank by sequentially traversing the element list, the last slave 
processor will possess the most ghost nodes as shown in right half of Fig. 2. In 
the upper left plot, elements with rank pi are traversed first and their corner 
vertices are assigned with rank pi marked with red squares. Elements with rank 
P 3 are processed at last. All the ghost nodes between CPU P 3 and p 2 and the 
ghost nodes between CPU P 3 and p\ will have rank P 3 . A special treatment is 
developed to improve the ghost nodes distribution as depicted in the lower right 
half of Fig. 2. We firstly stores the layer of ghost nodes into an array when they 
are first identified. Then alternates the ghost node rank once every other node 
to the secondly visiting processor’s rank. In the end, each processor will have 
nearly equal number of ghost nodes. 

This algorithm is developed solely for 3D surface meshes since most ghost 
nodes can only be shared by two processors. In volumetric meshes, ghost nodes 
can be on a curve which is shared by many partitions. Developing an algorithm 
evenly distributing the ghost nodes in volumetric meshes is in our future work 
list. 

4 Methodology 

4.1 Multigrid Explicit ODE Solver 

A novel time integration technique relying on multiple level meshes is used to 
compare with traditional single level explicit time integrator. Due to space lim- 
itation, we only describe the essence of the new algorithm. For further details 
please refer to [7,10]. This method is general enough to be embedded into dif- 
ferent modeling approaches besides FEM. The FEM formulation of our system 
can be found in [9]. 

Explicit time integration schemes have little computation cost per integra- 
tion step whereas small step size is required to accommodate system stability. 
By working on a single mesh, the schemes are limited to propagate spatial in- 
formation one layer at a step. Our proposed scheme targets this shortcoming 
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specifically and improve the propagation. Multiple meshes Mi’s can be gener- 
ated around the same geometry independently, where Mo has the highest res- 
olution. The number of elements N, in each Mi forms a geometric series with 
ratio q < 1. This is the key to have O(Nq) computation cost. There are three 
operators in MG, which either improve the solution at one level or transform a 
related problem onto another level. 

— The smooth operator G(-) is an explicit ODE solver, e.g. Forward Euler or 
Runge-Kutta family methods, who evolutes the state of Mi from time stamp 
t n to t n _ |_i- 

X^ 1 = G{X^,At) ( 1 ) 

— Restriction operator R(-) assigns the state of coarse mesh X„ ' with a linear 
combination of the state of the hosting fine element nodes in (2) . 



4°(p) 



R\_ iptf 1) ), \^LO^AX^\>e 

X%\p) + i2j_ 1 (AX’ (< - 1) ), | X)w (i-1) AX’ (i - 1) | < e 



where AX = X n — X n _\ and w’s are the weights. 

— Interpolation operator P(-), which is the inverse operation of R(-), interpo- 
lates Xn or AX ( ' l ' > to the next fine level in (3). 






PtHxP), \AX^(p)\>e 

xt 1 ] (p) + PT\AX®), < e 



Figure 3 shows the diagram of one V-cycle of MG time integrator. In contrast 
to single level iterative solvers, MG uses coarse grids to do divide-and-conquer 



Finest Level: 



Coarsest Level: 




Fig. 3. The diagram of one parallel MG V-cycle with 3 levels. 
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in two related senses. First, it obtains an initial solution from the fine grid by 
using a coarse mesh as an approximation. This is done by transforming the 
problem on the fine grid to the coarse grid and then improving it. The coarse 
grid is in turn approximated by a coarser grid, and so on recursively. The second 
way MG uses divide-and-conquer is in the frequency domain. This requires us 
to think of the error as a sum of eigenvectors of different frequencies. Then, 
intuitively, the work on a particular grid will attenuate the high frequency error 
by averaging the solution at each vertex with its neighbors. The nonlinear FEM 
engine we use has a linear computational cost of 300 FLOPS/node, for either 2D 
or 3D [9] . This is much greater than the cost of the restriction and interpolation 
processes in the multigrid algorithms, 2 x 42 FLOPS/node for tetralredra and 
2 x 30 FLOPS/node for triangle meshes. The overal computational cost of one 
MG V-cycle is dominated by C((Gi(-))) and yields 

G corn p = c ( N {i)) + Cost (mapping) 

< / m 0) ) + 84_/V (0) ]/(l - q) « 384JV 0 /(1 - «), 3D mesh 
" 1 [C(N { o)) + 60-/V (o) ]/(l - q) « 360-/V o /(l - q), 2D mesh 

The two sides of the above equation become equal as m approaches infinity. 
The increased computational cost is compensated by the stabler performance of 
multigrid integration, because in practice the multigrid algorithm can be inte- 
grated with a time step larger than twice that of the single level integrator. This 
is shown in the 2D examples in section 5. 

4.2 Parallel Multigrid Integrator 

The mesh at each level is partitioned into the same number of parts as in S^°b 
Each slave processor hosts a set of partitions, one of which is from each mesh 
level as shown in Fig. 4. In this example, a square mesh with N 2 nodes is 
considered. The next coarse level has N 2 / 2 nodes. Each slave possesses 1/8 of 
the mesh at each level. Ghost nodes information is exchanged. Let the cost of 
sending/receiving the state of one vertex as 7. Thus the communication cost at 
the finest level is 

C < ±m = K^+ l -N) 1 = A(V2+l)N 1 (4) 

Similarly, S*- 1 -* has the cost of 4(\/2+ \)N/\f2 r y. The overall communication cost 
Ccomm is bounded by 

/2 

Ccomm = ^2 (comm) < —j= — -[4(\/2 + l)iV]7 = (12v / 2 + 16)N^y (5) 

i v 

The communication dominates and eventually becomes a bottleneck when the 
slice size of coarse level is too small given the hardware and software dependent 
latency of preparing nodal message and the cost of sending such mesage. This is 
shown in Fig. 6. 
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Fig. 4. Each one of the 8 slave processors claims one slice of the mesh at each resolution 
level. 



4.3 Program Flow 

The main loop of the simulation consists of four steps: 

1. The front-end graphical workstation sends the latest user motion to the 
master node in the cluster via TCP/IP Socket. The master broadcasts this 
information to every slave. Each slave receives the new data. 

2. Each slave sends the current state of its ghost nodes to its adjacent slaves. 
Each slave also receives the data from its neighbors. Blocking receive is used. 

3. MG time integrator is invoked. The internal stress of the elements and the 
state of the vertices on each slave is updated and copied to a dedicated space 
in each slave’s local memory. 

4. Synchronizing among all the slaves ensures that every partition of the mesh 
has been updated before Pq collects the results from them and send the 
collection to the front-end client through TCP/IP socket. 

Notice that the simulation loop continues regardless of whether the master re- 
ceives new user motion from the front-end PC or not. In the last step, the 
simulation loop does not stop whether the front-end PC receives the new mesh 
or not. By this means, the cluster server has the capability to run the simulation 
independently. 

5 Results 

5.1 Comparing MG and Single Level Time Integrator 

MG integrator produces a more stable solution with little additional compu- 
tation cost, thus making our simulator viable for practical applications. The 



An Interactive Parallel Multigrid FEM Simulator 



131 




Fig. 5. Single level time integration versus. MG time integration. One node in a square 
mesh with 2628 elements with 8 partitions is lifted up. Left: The mesh is in the initial 
rest configuration with zero stress. Middle: Single level FEM update becomes unstable 
even under small deformation. Right: MG time integrator stays stable even under large 
user manipulation. 

parallelized version of MG time integrator has been implemented. The effect is 
demonstrated in Fig. 5. The single level explicit integrator did not propagate 
the stress concentration, due to user manipulation, fast enough into the rest 
portions of the mesh. The picture in the middle shows the mesh positions are 
diverging. The picture in the right was taken when using MG integration. The 
state of meshes at different resolution levels are updated simultaneously. MG 
time integrator stays stable even under large deformation. As indicated in (4), 
one MG cycle requires less than twice the cost of a single level method, so MG 
is cheaper to execute if it can stably simulate a dynamic system with the largest 
permissible time step max{hMG) > 2 max(h one _i eve i) . Because of the improved 
stability of MG, this is likely to be the case. 



5.2 Overall Performance 

To show the overall effect of parallelization, we ran experiments with a different 
number of processors with single level time integration, i.e. 2, 4, 8, 12, 16, 24, and 
32 slave nodes and different problem size, i.e. meshes with 321,668,1286, and 
2628 elements (Fig. 6). There are several observations that can be drawn from 
the experiments: 

1. For the same number of processors, r increases proportionally to the problem 
size. This verifies that the FEM engine [9] has linear cost. The slope of the 
linear increase is 0.069ms/element. 

2. For the same problem size, A spent by each slave node decreases as more 
slave CPUs are deployed. This verifies that the parallelization is efficient. 
Increasing cluster size does improve the update rate of the simulation. 

3. With the fixed number of processors, the master/slave communication cost 
increases as the problem size enlarges. For instance, using 12 slave nodes, 
each node spends 7ms to send data to the master node on the 321 elements 
mesh. On the 2628 element mesh, each node spends 32ms to talk to the root. 
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Fig. 6. The performance data is separated into groups. Each group associates with a 
different number of slave CPUs deployed during the test. Each group simulates the 
deformation of a square hyperelastic membrane under the influence of gravity with 
different discretization levels, namely, meshes with 321,668,1286, and 2628 elements 
are tested. The first bar in the left corner is labeled as 2(321). This means that there 
are 2 slave nodes used to simulate the deformation of the mesh with 321 elements. 



6 Conclusion 

The proposed system serves as the prototype of a scalable distributed computing 
environment for interactive deformable object simulation applications. Explicit 
MG time integrator is parallelized as well to obtain stabler performance com- 
pared to single level explicit integration schemes. 

The current system implemented the parallelization of FE modeling and 
Multigrid time integrator on 3D surface object. Our parallelization can be ap- 
plied to 3D volumetric meshes by just replacing triangular element stress formula 
by tetrahedron element stressformula. The element assembly, load balancing, 
MG integration are the same. Recall that in section 3.1 the ghost node distribu- 
tion algorithm is not optimized for volumetric meshes. Future investigation will 
develop algorithms for evenly distributing ghost nodes in volumetric objects. 

Haptic interface was not included in our cluster system though not limited 
to. The main issue is the non-static frame rate on the graphical front-end client, 
which is caused by public TCP/IP connection between the cluster server and 
the graphical client. In order to incoporate haptic feedback, the frame rate fluc- 
tuation need to be accurately predicted. 

This system is cost effective since it is built from Linux PCs and commercially 
available hardware. The system is also highly scalable, because replacing an 
old or adding a new computing node only effects the cluster locally. The load 
balance is well achieved by METIS off line. MPI and C/C+- 1- on Linux are 
our only programming environment which is portable. Additionally, the data 
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distribution is independent with the physical network topology since network 
devices themselves handle the relay of packages between nodes are not directly 
connected. 

The inter-processor communication pattern is very complex in the heteroge- 
neous cluster system. The communication delay is not uniform across the cluster. 
The impact of time delay on user/haptic interaction will be studied in the future. 
Further performance experiments and data analysis need to be carried out in 
the future in order to locate the shortcomings and bottlenecks in the system. 

The purpose of the simulator is to introduce a successful prototype of a 
scalable interactive distributed computing system used for deformable object 
simulation applications. The proposed system design, partitioning schemes, and 
the interactive client/server architecture will bring further research into tele- 
surgical simulation, remote scientific education, and interactive large scale data 
set visualization. 
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Abstract. The Extended Finite Element Method (XFEM) is a tech- 
nique used in fracture mechanics to predict how objects deform as cracks 
form and propagate through them. Here, we propose the use of XFEM 
to model the deformations resulting from cutting through organ tissues. 
We show that XFEM has the potential for being the technique of choice 
for modelling tissue retraction and resection during surgery. Candidates 
applications are surgical simulators and image-guided surgery. A key 
feature of XFEM is that material discontinuities through FEM meshes 
can be handled without mesh adaptation or remeshing, as would be re- 
quired in regular FEM. As a preliminary illustration, we show the result 
of XFEM calculation for a simple 2D shape in which a linear cut was 
made. 



1 Introduction 

Neurosurgeons plan surgery from patients’ structural and functional images. 
During surgery, neuronavigation systems display the positions of surgical ins- 
truments in preoperative images. However, the brain deforms in the course of 
a surgery. These deformations occurs principally following the opening of the 
dura, the drainage of the cerebrospinal fluid (CSF), the retraction of tissues and 
the successive resections of, say, a tumor [5]. The usefulness of neuronavigation 
systems is then limited: the current brain shape no longer corresponds to that of 
the preoperative images. Intraoperative image acquisition can partially circum- 
vent this limitation by capturing the new configuration of the brain, but such 
image acquisition is limited in signal-to-noise ratio and spatial resolution by the 
time constraints of the surgical procedure. Additionally, not all imaging moda- 
lities (particularly functional ones) are available intraoperatively. Consequently, 
intraoperative monitoring and surgical navigation can be significantly improved 
by estimating the deformation of the brain and projecting preoperative imaging 
data into alignment with the subject brain. 



S. Cotin and D. Metaxas (Eds.): ISMS 2004, LNCS 3078, pp. 134—143, 2004. 
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Nonrigid registration techniques are numerous. One approach is to use biome- 
chanical models to encapsulate the mechanical properties and behavior of the 
brain. We use intraoperative MRI images in order to compute the displacements 
of the cortical and ventricular surfaces of the brain. The displacement field then 
drives the brain model in place of the forces. Deformations throughout the brain 
are calculated using the finite element method (FEM) with a linear elastic be- 
havior law [18]. Our approach and implementation are similar to those of Ferrant 
[5] [7] [6], 

Virtually all past studies on brain deformation modelling have focused on 
the brain shift at an early stage of the procedure, before significant resection 
has taken place. The precision achieved in the prediction of displacements in 
this particular context are good, with for instance, landmark matching errors 
of 0.8 ± 0.4 mm at the surface of the brain and 1.1 ± 0.7 mm for the interior 
( me an± standard deviation) [7], which is comparable to the size of the voxels. 

In comparison, modelling of retraction and resection is still in its infancy. 
There have been limited investigations of the forces involved in these two sur- 
gical tasks [10]. There has also been an effort to create a so-called “smart re- 
tractor” capable of measuring forces intraoperatively; this device has already 
shown promising results [1] . Intraoperative images will be very useful when used 
in conjunction with this device. Ferrant et al. [6] has used intraoperative MRI 
and captured brain deformation in the presence of resection, but the model of 
resection simply consisted of clipping the deformed brain with the resection cav- 
ity. 

Conventional FEM has serious limitations for modelling the tissue disconti- 
nuities associated with the resection of a tumor. The displacement field that 
is the solution of the finite-element (FE) calculation must be continuous inside 
each FE. Furthermore, the displacement at any node may take on only one 
value. Consequently, the modelling of a discontinuity requires that discontinuity 
boundaries be aligned with boundaries of the elements and that nodes lying on 
these boundaries be duplicated. 

Extensive work has been performed in the domain of surgical simulation on 
the problem of cutting of a FE mesh [13] [3]. Most of the proposed solutions and 
algorithms are based on a subdivision method [2] [9] [12] [8] [13]. All elements inter- 
sected by the cut are divided into sub-elements in order to create a boundary of 
finite elements aligned with the cut. The subdivision is subject to the constraint 
of a good aspect ratio for the new elements. 

The main drawback of this method is the rapid growth of the numbers of 
nodes and elements in the mesh. In addition to the subdivision calculation, the 
larger mesh size increases the computation time, and it is challenging to maintain 
computationally efficient parallel data structures as the mesh evolves. Of course, 
it is important to keep in mind that real-time performance is absolutely essential 
for surgical navigation and simulation. 

Alternative methods have been evaluated to avoid these dramatic changes in 
the mesh. For example, the mesh may be adapted to the geometry of the cut: 
some nodes are selected, and they are then relocated to cling as best as possible 
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to the cut geometry. However, an offset can remain between the boundary formed 
by these relocated nodes and the cut. Element degeneracies can also happen [14]. 
Depending on the method used, the distortion of the mesh can produce elements 
with unacceptably large aspect ratios. One solution is a remeshing, but this will 
in turn lead to an increase in the computation time [15]. 

We propose a new method for cutting meshes in arbitrary ways without mesh 
adaptation or remeslring, thereby avoiding all of the above drawbacks. This 
method allows the object to be modeled by finite elements without explicitly 
meshing the cut surfaces. Discontinuities can then be arbitrary located with res- 
pect to the underlying FE mesh. In addition, no remeslring is required when the 
discontinuity changes shape. Other appealing features of the method are that 
the FE framework and its advantages (sparsity and symmetry) are retained and 
that a single-field (displacement) variational principle is used. This method is 
called “extended finite element method (XFEM)” and was introduced in 1999 in 
the field of fracture mechanics for the study of crack and related failures, such 
as for bridges and airplanes [11]. 

The paper is organized as follows. In Section 2, we review the basic principles 
of FEM. In Section 3, we discuss the theory of XFEM. In Section 4, we provide 
a proof-of-concept simulation for a simple 2D object. Finally, in Section 5, we 
discuss the role and benefits of XFEM for handling the deformation of the brain 
in surgery, especially in the presence of significant retraction and resection. 



2 Review of Basic FEM Principles 

The problem of finding the displacement field u(x) such that the weak form 
of the equations of linear elasticity is satisfied is equivalent to determining the 
displacement field which minimizes the total deformation energy E 

E=— I u e cLQ — J bucLQ— j tudr. (1) 

2 J o Jo J r t 

The quantities in (1) are as follows. e(x) and cr(x) are the strain tensor and the 
stress tensor, respectively. b(x) is the body force applied to the solid, while t(x ) 
is the traction force applied to its surface. Q represent the volume of the solid 
and r t represents the surface of the solid on which traction is applied. 

To solve the linear elastic problem, we need to discretize the equations. In 
particular, we need an approximation u h for the displacement field u. The FEM 
approximation is defined by 



N 

u h (x) ='^2ip i (x)u i , (2) 

i-1 

where the uj s are the discrete unknowns to be determined and the ipj s are 
basis functions, called shape functions. These functions must obey 2 conditions. 
First, tpi has a compact support oji, which corresponds to the union of element 
subdomains connected to node i. Second, we have 
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1 if i = j 
0 if i^j 



on u>i , 



where the Xi’s, i = 1, N, are the coordinates of the nodes. 



( 3 ) 



Equations (2) and (3) yield the following property 



N 

u h (xi) =y^<p i (xi)u i = u-i. (4) 

i=l 



The FEM unknown iq can be shown to be the displacement field value at the 
node x-i. The FEM displacement field interpolates nodal displacements. 

Finally, the introduction of the FE approximation (2) in the minimization 
of (1) leads to the following system of linear equations 



where 



Ku = f 



FT; II 



= fi i,j = 1, 
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( 5 ) 



( 6 ) 



( 7 ) 



and H is Hooke’s tensor. The final step is to solve (5) for the displacement Ui. 
One can then use (2) to align preoperative and intraoperative images. 



3 Introduction to Basic XFEM Principles 

3.1 Fundamental Equations 

This section introduces the fundamental ideas of XFEM. For details regarding 
the theory and various implementation issues, the reader should consult, e.g., 
[11] [17] [4] [16], 

The key of this method is to create a new displacement-field approximation 
by enriching the FE approximation (2), that is by multiplying some of the FE 
nodal shape functions by discontinuous functions. This enrichment can be made 
to take local form by only enriching those nodes whose support intersect a region 
of interest. We have 

n Ei 

u h (x) = y; <pj(x)uj + y: ipi(x) y^g j (x)a ji . 

iel iej j — 1 



(8) 
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The quantities in (8) are as follow. The ipf s are the FE shape functions and 
the gj' s are the XFEM enrichment functions. We denote by I the set of all N 
nodes in the domain, and by J the subset of / corresponding to the n E enriched 
nodes. Ui and dji are nodal DOFs and n El denotes the number of enrichment 
functions for node i. The additional DOFs are associated with nodes that 
are enriched. 

An important consequence of the XFEM function enrichment is that the 
approximation does not interpolate nodal displacements for enriched nodes Xi, 
i.e., 

n Ei n Ei 

u h (xi) = Y,'Pi(x i )u i + Y i 'Pi{xi)Y,giixi)a ii = Ui + ^2gj(xi)aji u t . (9) 

iel iej j — 1 j = 1 



3.2 Choice of Enrichment Functions and Enriched Nodes 



We denote by Fj the crack surface. Any function that is discontinuous across r r j 
can be used to model an arbitrary discontinuity in u h (x) . The simplest choice is a 
piecewise-constant function that changes sign at the boundary F ( i, the Heaviside 
function 



( 1 for ( x - x*).e n > 0 
\ -1 for (x - x*).e„ < 0 



(10) 



where a; is a sample point of the solid, x* is the point on the crack that is the 
closest to x , and e n is the outward normal to the crack at x* 1 (Fig. 1(a)). 
The nodes that are enriched by this function are those for which the support 
intersects the crack. 




nuity (2D case), (b) Local coordinates of the crack-tip enrichment functions (2D case). 

However, this function is not sufficient to model accurately the tip of the 
crack when this tip terminates inside an element. Indeed, the node and so all 
of its support is enriched. Consequently, the crack would be modeled as though 
the “crack-tip” was extended till it intersects the element edge. 



1 



Outward” is defined in an obvious way based upon the relative positions of x and I'd- 
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Consequently, nodes whose supports containing a crack-tip are not enriched 
with the Heaviside function, but with specific crack-tip enrichment functions that 
ensure that the crack terminates precisely at the location of the crack-tip. The 
crack-tip enrichment relies on functions that incorporate the radial and angular 
behavior of the asymptotic crack-tip displacement field, which is two-dimensional 
by nature (Fig. 1(b)) 

{Fi(r,9)}t = i = {y/fsin(^),Vfcos(^),Vfsin(^)sin(Q),y/rcos(^)sin(6)}, ( 11 ) 

where r and 6 are the local polar-coordinates 2 . 

The XFEM approximation for a single pair of crack and crack-tip is thus 

N 4 

uh ( x ) = Y Vi(x)ui + y <pj(x)H(x) aj + y M x )CY F ^ x ) c k)- ( 12 ) 

i—1 jej keK l—l 

The quantities in (12) are as follows. The Uj’s are the nodal DOFs associated 

with the continuous part of the FE solution, the cij's are the nodal enriched 
DOFs associated with the Heaviside function, and the c[.’s are the nodal enriched 
DOFs associated with the crack-tip functions. I is the set of all nodes in the mesh. 
J is the set of nodes whose shape function support is cut by the crack interior. 

K is the set of nodes whose shape function support is cut by the crack-tip (x c ). 

With D denoting the crack geometry, we thus have the formal definitions 

K = {k e I : x c e Wkj J = {j e I : u>j ft D yf 0, j ^ A'}, (13) 

where Uk denotes the compact support of the node k and Wk its closure. The 
above equations can easily be generalized to several pairs of cracks and crack- 
tips. 

To obtain the discrete XFEM equations equivalent to the FEM equations 
(5)-(7), we must substitute the approximation expression (12) in the total-energy 
expression (1) and minimize the resulting expression. 

Additional details regarding the equations can be found in [17]. 

While FEM requires a remeshing and the duplication of nodes along the 
crack to take into account any discontinuity, the XFEM requires identification 
of nodes belonging to the sets J and K and the computation of the stiffness 
matrice with enrichment fonctions. Because of added nodal DOFs in XFEM, 
stiffness matrix is larger than in FEM. 



4 Proof-of-Concept 2D XFEM Simulation 

To evaluate the abilities and potential of XFEM for surgical guidance and simu- 
lation, we have performed preliminary tests on simple 2D objects such as rectan- 
gles and ellipses containing a line-segment crack discontinuity. An exploratory 

2 The first function is discontinuous on the crack faces. 
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Fig. 2. (a) Mesh and crack geometries before deformation, (b) Results of XFEM when 
displacements (14) and (15) given in text are applied along the crack. 



program was written in Matlab. The Matlab PDE toolbox was used, but only 
to initialize a triangular mesh from a domain boundary. 

The inputs to the program are the mesh definition and the crack geometry 
(Fig. 2(a)). One begins by identifying the mesh elements that are fully intersected 
by the crack and the mesh elements that contain a crack-tip. The number of 
DOFs for each node is then defined: 2 for a non-enriched node, 4 for a Heaviside- 
function-enriched node or 10 for a near-tip-functions-enriched node. As with 
FEM, each elementary stiffness matrix is computed, taking into account the 
enrichment functions, and the global stiffness matrix is subsequently assembled. 
The application of force or displacement constraints is performed in similar ways 
in both FEM and XFEM. 

To illustrate the method, we have chosen to constrain the displacement 
along the (linear) crack. We have applied a displacement of (0.0266, 0.0103) and 
(—0.0266, — 0.0103) 3 to each intersection, with the mesh, of the right and left 
crack lips, respectively. The intersection of the crack with the element containing 
the crack tip was left free (i.e., no displacement was imposed) to avoid having 
excessive constraints near the tip. Indeed, we have noticed that element flipping 
can sometimes happen in this situation. 

The application of this displacement constraint is straightforward in XFEM. 
For the intersection (* mt ) of an element defined by 3 nodes enriched with the 
Heaviside fonction, the relation between the nodal DOFs 4 are 

f 1^x1 + O'xl T ^2 x2 T t p'2(p'int) ^x2 = 0.0266 /-..\ 

( T ^Pl ip^int) & yl T '^y2 T ^f2 ®y2 — 0.0103, 



3 Each point is of the form (x,y), with the x and y axes being respectively horizontal 
and vertical in Fig. 2. 

4 We consider the intersection lying on the element boundary between node 1 with 
DOFs (u x i,Uyi,a x i,a y i) and node 2 with DOFs (u X 2,u y 2,a X 2,a y 2). 
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for the right lip where the Heaviside function (10) is equal to 1. Similarly, we 
have 

f Pi Hxl Pi i&int) &xl T P2{p^int) U x2 P2 (&int) &x2 — 0.0266 / , r \ 

\ Plip^int) 'Uyl Plip^int) & yl T P2 {p^int) Uy2 P2 ip^int) & y2 — 0.0103, 

for the left lip where the Heaviside function (10) is equal to —1. The result of 
these displacement constraints is an opening of the crack, which is illustrated in 
Fig. 2(b). 

This example confirms that XFEM can elegantly and efficiently take into ac- 
count (crack) discontinuities in the study of the mechanical properties of objects. 
In particular, remember that no mesh adaptation or remeshing is required. In 
contrast with FEM, solution displacement fields can now contain discontinuities 
inside the finite elements. Observe that the triangles that appear to have been 
added in Fig. 2(b) were added for display purposes only. They are needed to 
show the new boundary of the crack. However, none of these additional nodes 
or elements are involved in the XFEM calculation. 

5 Conclusions and Perspectives 

XFEM is particularly well adapted to deal with the general problem of cutting 
through a 2D or 3D finite-element mesh. This is required to deal with disconti- 
nuities, i.e., cracks or cuts. The main feature of XFEM is that it can deal with 
discontinuities without having to perform computationally-expensive mesh adap- 
tation or remeshing. Note that the technique applies to multiple discontinuities 
that have arbitrary locations and shapes. 

The implication for surgical guidance and simulation are clear and signifi- 
cant. In surgery, XFEM could be very useful in the modelling of retraction and 
resection, each of these surgical procedures inducing discontinuities in tissues. 
The pieces of information we need are the discontinuity geometry and the dis- 
placement contraints along the organ surfaces and the discontinuity boundary. 
The incision surface allowing the insertion of the retractor can be determined 
by tracking. It can also be inferred from an intraoperative MRI image show- 
ing the retraction pathway. The displacements caused by the retractor can be 
calculated from distances between segmented brain boundaries along the retrac- 
tion path and the calculated incision surface. The modelling of resection is more 
complicated given that the brain can swell during this surgical procedure and 
that this swelling is not visible in intraoperative images. However, the difficult 
task of removing finite elements according to the boundary of resected areas 
can be done accurately with XFEM. Indeed, all elements falling entirely in the 
resected area can be removed. With the remaining elements, we can precisely 
specify the boundary of the resected cavity by adding discontinuous functions 
to nodes located along this boundary, which allows us to cancel the presence of 
the elements on the resected side of the boundary. 

The next step in our XFEM work will be the application and evaluation of 
the ideas and techniques proposed above to MRI images of a phantom submitted 
to retraction and resection. 
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Abstract. We have earlier proposed a voxel-based representation of an elastic 
object that can respond to a user’s input in real-time for haptic rendering. We 
called it shape-retaining chain-linked model or S-chain model. The S-chain 
model is constmcted from the 3D voxel data of an object, where each voxel is a 
chain element that is linked to its six nearest neighbors. Its deformed configura- 
tion is computed, upon the user’s input, by propagating outward the unabsorbed 
input force from the interaction point as if a chain is pulled or pushed. The de- 
formed configuration is then used to compute its disturbed internal energy that 
is reflected to the user. The basic idea of force rendering is that the reflected 
force is proportional to the number of chain elements that are displaced from its 
initial position. This simple nature of the model allows very fast deformation of 
a volumetric object so that it can be utilized in real-time applications. In this 
paper, we present a mechanical interpretation of the haptic model as to how the 
reflected forces are being computed by utilizing spring-and-cylinder units. Fur- 
thermore, we investigate the quality of the haptic feeling of the S-chain model 
by comparing with that of FEM against the human haptic perception. The result 
of our experiments demonstrates that S-chain model provides not only the real- 
time performance but also the quality of FEM with respect to our haptic sense. 



1 Introduction 

The utility of surgery simulation in medical education is increasing with advances of 
the technology in virtual reality. Especially with a new trend of minimally invasive 
surgery (MIS), where the surgical tools are inserted into small poles and surgeons 
watch the surgical area displayed on a monitor via a scope, the need for training tools 
is growing fast. At the same time, since the MIS restricts the surgical environment, it 
becomes more feasible to develop the matching virtual environment for the simula- 
tion. For an MIS surgery simulation, we must have a virtual object representing the 
bodily organ or some portion of the body, which interacts with a user mostly via tools, 
and realistic feedback from the virtual object upon interaction. The virtual object is 
usually deformable and desirably volumetric. To have a realistic feedback, the sen- 
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sory mode should include not only the visual but also the haptic which provides the 
kinesthetic information. Moreover, the response should be in real-time: within 30Hz 
for visual feedback, and within 1000Hz for haptic feedback. 

Previously, we have proposed a voxel-based representation of an elastic object that 
can respond to a user’s input in real-time for haptic rendering [1-4]. We called it 
shape-retaining chain-linked model or S-chain model. The S-chain model is con- 
structed from the 3D voxel data of an object, where each voxel is a chain element that 
is linked to its six nearest neighbors. Its deformed configuration is computed, upon 
the user’s input, by propagating outward the unabsorbed input force from the interac- 
tion point as if a chain is pulled or pushed. The deformed configuration is then used to 
compute its disturbed internal energy that is reflected to the user. The basic idea of 
computing the reflected force is that the reflected force is proportional to the number 
of chain elements that are displaced from its initial position. This simple nature of the 
model allows very fast deformation of a volumetric object so that it can be utilized in 
real-time applications. 

In this paper, we present a mechanical interpretation of our haptic model as to how 
the reflected forces are being computed by utilizing spring-and-cylinder units. In our 
early work of [1], we incorporated force- voltage analogy (duality) concepts in order 
to develop a haptic model from the S-chain construction. As a result, the serially- 
connected capacitor model was introduced. However, since there were some doubtful 
comments on the use of ‘capacitors’ in the haptic representation, we came up with 
another interpretation that is more appropriate in mechanical sense, using spring-and- 
cylinder units. We would like to explore the spring-and-cylinder haptic interpretation 
for our S-chain model. Furthermore, we have investigated the quality of the haptic 
feeling of the S-chain model by comparing with that of FEM against the human haptic 
perception. The result of our experiments demonstrates that S-chain model provides 
not only the real-time performance but also the quality of FEM with respect to our 
haptic sense. 



2 Shape-Retain Chain Linked Model (S-Chain Model) 

Without loss of generality, we will use ID case for an easier discussion of the model. 
The 2D and 3D cases are the natural extension of the ID case as described in our 
previous work [1]. The topmost configuration in Figure 1(a) shows the initial configu- 
ration of the ID model having 5 nodes. Each node is connected to its nearest 
neighbors, and it has a finite range of free motions without disturbing its neighbor. 
The range is drawn like a chain element for each node as shown in the figure. (Note 
that the range shown in Figure 1 is 2D rather than ID, but its degree of freedom is 
only ID moving to the left or to the right.) Suppose a user selects the rightmost node 
(or e chain element) and pulls it to the right. The e will be moved to its right up to its 
maximum stretch limit as shown in the 2 nd configuration of Figure 1(a). If the user 
continues to pull e to the right, e 2 is moving together with e upon reaching the 
maximum stretch limit for e as shown in the 3 rd configuration of the figure, and e 3 
will be moving together with e 2 and e, upon reaching the maximum stretch limit for e 2 , 
and so on. These are the basic steps of computing the deformed configuration of S- 
chain model. Given an input force, the force is absorbed by each chain element dis- 
placed. If there remains an input force, the linked chain elements will absorb the force 



146 



Jinah Park, Sang-Youn Kim, and Dong-Soo Kwon 



e, e 4 


e ? 6, 


e 5 e 4 


“3 




1 ^ 


e< e„ 


3^ (at 




1 1 * 


e 5 e 4 


e i e 7 e. 


1 


1 1 * 




Xinit|_^k4 


<k, <k. 


J 

e< e„ e 


^ S x > 

f ^ma) r 

1 (tit (b) 

[ e 2 [ | e l ] 


GY] 


| c) (c) 




1=1 



Fig. 1 . ID S-chain Model having 5 nodes 

as much as possible by displacement. This propagation continues until there is no 
more input force to be absorbed. Since the process is one-way out and each computa- 
tion is no more than simple comparisons, the computation is very fast. Also, since the 
computation is done locally, the algorithm does not depend on the size of input (the 
number of nodes comprising the object). 

The significant difference between 3D ChainMail algorithm (originally proposed 
by Gibson [5]) with our S-chain model [ 1 ] is that our algorithm computes the defor- 
mation against the initial resting state where the internal energy is in equilibrium. This 
condition forces to retain the deformed shape under the same force; and, therefore, 
allows us to relate the perturbed internal energy of the deformed object to a haptic 
feedback force directly and solely based on its deformed configuration. 
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Fig. 2. 3D S-chain Model having 75x75x75 nodes 

Figure 2 shows a deformed configuration of a cubic object having 75x75x75 
nodes. We display here only the connections of nodes. An arrow is overlaid to show 
the interaction path and the white sphere indicates the interaction point with the user 
(via a haptic device). Although only the surface elements are displayed in Figure 2, 
the model is truly 3D and the deformation computation is done in 3D in real-time. 

The deformed configuration is then used to compute its disturbed internal energy 
that is reflected to the user as a haptic feedback. The haptic rendering is performed 
based on the simple observation that the reflected force is proportional to the number 
of chain elements that are displaced from its initial position. Each chain element has 
its own ‘material property’ that is assigned according to original voxel classification. 
This stiffness parameter k, and the ‘appropriate’ displacement id) of the chain element 
are then used to compute the force using Hooke’s law: F. = k d. The force reflected to 
the user is the summation of the force for all chain elements displaced: F = X F.. The 
next section describes what we mean by the appropriate displacement, and how they 
are transferred to the mechanical representation. 



3 S-Chain Model: Mechanical Representation 

Figure 1(b) shows the internal energy of the ID model in Figure 1(a) at the state 
where three chain elements (e r e 2 and e 3 ) are displaced and where two (c and e 2 ) of 
which are stretched to their maximum. The chain element e 4 is drawn together to 
show the reference of the initial position. A spring is used to capture the potential 
energy absorbed by each chain element. The magnitude of the reflected force to be 
generated in this sample case is the sum of the individual forces generated by e p e 2 
and e F = k, X„„ + k 2 X„, + k, X H where k.’s are the stiffness coefficients for each 
node (or chain element), X max is the stretch displacement limit for the node, and X d is 
the displacement of the 3 rd node e 3 . Assuming that X max is the same for all nodes in this 
case, we can generalize the force computation as follows: 

F = (i>,X ma J + k n x„ (!) 

i=l 

As shown in Figure 1(b), the springs are connected in parallel to represent the system. 
For further discussion on why the springs should be connected in parallel, not in se- 




148 Jinah Park, Sang-Youn Kim, and Dong-Soo Kwon 



ries, the reader is referred to [1], And each spring has its maximum stretching limit 
(and compression limit), regardless of the absolute displacement of the node positions 
from its initial position. Only those nodes that are displaced will contribute to the 
reflected force. In order to capture all of these, we present a mechanical unit having a 
spring and a frictionless cylinder, connected as shown in Figure 1(c). The frictionless 
cylinder acts as a switch to activate the spring. The length of a cylinder indicates the 
buffer length as to how much it awaits to activate the spring connected to it. Let the 
node that is in direct contact with the user (or a haptic device) be the USCE (user 
selected chain element). Then the cylinder connected to the USCE has a zero length to 
indicate that the spring is activated immediately. The spring-cylinder units are con- 
nected in parallel from the USCE to the rest of nodes in the model. In Figure 1(c), 
four of such units are connected from e, to e 2 , e 3 , e 4 and e 5 . 

The parallel connection is simple to understand in ID case as shown in Figure 1(c). 
Let us now consider 2D case (of 8x6) as shown in Figure 3. The USCE (displayed as 
a big circle) is the node that is in direct contact with the user. Let us define the ‘major 
linked-chain’ as those chains that are directly linked to USCE. They are indicated in 
thick lines in Figure 3. Those nodes that are on the major linked-chain are called 
‘leading chain elements’ in our discussion. Note that there are three major linked- 
chains in the case shown in Figure 3, and if the USCE is somewhere in the middle, 
there will be four major linked-chains. 

The parallel connection of spring-cylinder units is shown in Figure 4, where each 
line indicates an individual spring-cylinder unit. The strategy of connecting the unit is 
first to connect the leading chain elements directly to the USCE; and then to connect 
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Fig. 3. 2D S-chain Model 




Fig. 4. 2D S-chain Model of Figure 3 with spring-cylinder connections 
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each leading chain element to its leading chain elements. The dark lines in Figure 4 
are the links directly connected to the USCE, and the light lines are the links con- 
nected to the leading chain elements. The strategy is the same for the 3D model, 
where the leading chains are extended to the 3D space (into 6 directions in most gen- 
eral case). 

This (parallel) connection is represented by a simple summation in equation (1). 
Out of all the spring-cylinder unit connections, only those springs that are activated 
by the cylinder of its pair contribute to the force sum-up. 



4 Haptic Output Quality Evaluation 

In order to evaluate our haptic rendering algorithm, we experimented with a 3D liver 
model constructed from CT data. For the experiments, we used a PHANToM™ de- 
vice and an 800MHz Pentium III processor with 512MB of main memory. The soft- 
ware program is coded in VC++ with OpenGL, and the simulation is performed at the 
haptic update rate of 1000Hz and the graphic update rate of 100Hz. Figure 5 shows a 
picture performing our experiment. Figure 6(a) and (b) show graphical representations 
of the deformed configuration of the liver at the compressed state and the stretched 
state, respectively. 




Fig. 5. Experiment 




Fig. 6. 3D chain model of a liver 
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To investigate the haptic behavior of our S-chain model, we conducted a compara- 
tive study. We plotted the required force related to the user’s input when the object is 
modeled with the S-chain model and with FEM. According to the basic science re- 
search of [6], the behavior of a human liver is non-linear and anisotropic, and the non- 
linear shape of the stress and strain curve can be well approximated by a 3 rd degree 
polynomial curve as in Equation (2) where a is a compressive stress, 8 is a compres- 
sive strain, a t and a 2 are non-linear factors, and E is Young’s modulus. We set the 
stiffness value of the S-chain model and the FEM as non-linear following the equa- 
tion. 

C = £8 (1 + £q8 + fl 2 8 2 ) (2) 



7-i 




7-i 


6- 


j 


6 


§5 


J 


§5 




<D 




a» 




FEM 




LL. 


/ 


U- 


-o 3 




TJ 3 


a> 


■ 


<D 


t> o 

<D Z 


¥ 

J 


o 2 
a> L 


5= 






a> 1 

a 1 




a* a 

a: 1 


0 




0 




Amount of deform<ition(x1000) 



1 2 3 4 5 

Amount of Deform ationfxl 000) 

(a) (b) 

Fig. 7. Haptic Force Outputs from FEM and S-chain Model 



Figure 7 (a) and (b) show the haptic results when the human liver is modeled with 
FEM and the S-chain model, respectively. The amount of deformation refers to the 
number of voxels that were displaced. In the FEM simulation, FEM model with 
commercial program Abaqus was used. As observed in two graphs of Figure 7, the 
characteristics of the curves are very similar to each other. It took approximately 3 
minutes of computation time to obtain the reflected force per unit sampling step with 
the Abaqus 1 , while it took less than 0.01 seconds with S-chain model. Each dot on the 
graph shows each sampling step. 

The differential threshold for the force that the human can reliably discriminate is 
about 10% over a force range of 0.5-200N and the threshold increases to 25% forces 
smaller than 0.5N [11]. We calculated the discrepancy error of the S-chain model 
against the values obtained from the experiment with FEM, and verified that the dif- 
ference between the reflected force from the S-chain model and that from the FEM is 
smaller than the differential threshold that the human can reliably discriminate. 



1 Several real-time deformable modeling methods with FEM have been proposed in the past, 
such as in [7-10], However, we were not able to perform the comparable study with the algo- 

rithms. 
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Fig. 8. Discrepancy of values computed from the two graphs in Figure 7 



5 Summary 

In this paper, we present a mechanical interpretation of the haptic model for S-chain 
construction as to how the reflected forces are being computed by utilizing spring- 
and-cylinder units. Furthermore, we investigate the quality of the haptic feeling of the 
S-chain model by comparing it with that of FEM against the human haptic perception. 
The result of our experiments demonstrates that the S-chain model provides not only 
the real-time performance but also the quality of FEM with respect to our haptic 
sense. In medical simulation, the use of volumetric data, which are readily available 
from modern medical imaging scanners, can bring much advantage to the community. 
We believe that our proposed S-chain model as an alternative modeling scheme may 
build a new venue for advancement in virtual reality in medicine. 
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Abstract. The paper describes three new tissue deformation 
algorithms. We present a Mass-spring simulation with a quasi-static mod- 
ification of the Euler integration to increase the stability of the simula- 
tion. A directed length correction for springs and an algorithm called 
Dragnet are suggested to enhance propagation of large local displace- 
ments through the Mass-spring mesh. The new algorithms are compared 
with methods already in use. The combination of Dragnet and the quasi- 
static Mass-spring modification is used for the interactive real-time sim- 
ulation of an ophthalmological procedure, the removal of the Internal 
Limiting Membrane (ILM). 



Introduction 

The Internal Limiting Membrane (ILM) is part of the retina of the human eye. 
Under certain circumstances, a surgeon must remove the ILM to get access 
to the lower parts of the retina. The removal of the ILM is a difficult task 
because retina damages must be avoided to save the eyesight of the patient. A 
training module was developed to train ILM removal. Platform for the module 
is the commercial ophthalmosurgical simulator EyeSi 1 )!]. Different simulation 
algorithms are compared to find a suitable approach for modelling the rigid 
tissue of the Internal Limiting Membrane. 

Previous Work 

The Finite Element Method (FEM) [2] [3] and Mass-spring models [2] [3] [4] [2] with 
explicit or implicit integration [5] are common to simulate soft tissues. FEM and 
Mass-spring with implicit integration are often used to model rigid tissue because 
of their numerical stability. Both approaches require inverting a large sparse ma- 
trix as pre-processing step. After topological changes like cutting and tearing or 
modifications of the fix-point constraints the expansive pre-processing step must 
be carried out again. For interactive real time applications considering cutting 
and tearing, Mass-spring with explicit Euler integration is still a valid approach. 
The algorithm is fast, so that the mesh quickly reacts to interactions. Topologi- 
cal changes like cutting and tearing do not need additional calculation time for 

1 EyeSi is a product of the VRmagic GmbH, Mannheim, Germany. www.VRmagic.com 
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changes in the simulation structure. The algorithm is simple and easy to imple- 
ment. On the other hand, propagation of large local displacements is a weakness 
of that method since the translation of a node only affects its direct neigh- 
bours within one time-step. Another problem is the simulation of rigid tissue. 
It requires short time-steps for numerical stability that increases the calculation 
time. 

The following equations describe the explicit Euler integration step for a node 
in a mass-spring simulation. 






Fj(t){x,v) Af 

mi 



(1) 



Ax l{ t) = V i(t y At = (u i(t _ 1} + Avi( t) ) • At (2) 

where ( t ) is the time-step, m,; is the mass and x , v and F^ t \(x,v), are the 
position-, the velocity- and the force- vector of node i. 

The Force F^(x, v) usually has the form 

Fi{t)( x i v ) = Ri(t)( v ) + + Ei (3) 

where R^ t \(v) is the friction force, S^u^x) is the sum of the spring forces and 
is an external force to the node. 

Many modifications of Mass-spring exists to increases the stability and the 
propagation of large local displacements. 

Provot[6] suggests a correction of spring lengths if they exceed a critical 
maximum or minimum value. In a pre-processing step of Mass-spring the nodes 
of every spring exceeding its limits are moved along the spring direction to 
keep the spring in range again. This accelerates the propagation of large local 
displacements through the mesh. The length correction distributes local stress 
and increases numerical stability. Provot uses the length-correction to model 
rigid cloths behaviour. 

The standard Euler algorithm calculates forces for every node in the mesh 
and then the displacements. Brown et. al.[7] suggest to calculate the node-force 
and displacement node by node. To fasten propagation the order in which the 
node processing is done depends on the interaction point. The closer the nodes 
are to the interaction point the earlier they are processed. Brown et. al. use the 
algorithm to calculate the deformation of blood vessels. 

Mosegard[8] combines the two approaches of Brown and Provot to the LR- 
Spring Mass model. He uses LR.-Spring Mass for a cardiac surgical simulation. 

Chain-mail proposed by Gibson [9] is an algorithm that is well suited for real- 
time simulation of plastic deformations of volumetric bodies. Instead of using 
forces and elasticity Chain-mail handles the displacements of the nodes directly. 
In a relaxed mesh, every node has a certain range in which he can move with- 
out affecting its neighbours. If the distance to one of his neighbours exceeds a 
minimum or maximum limit, the node shifts its neighbour with him to keep a 
valid distance. The Chain-mail method calculates the movement for the x-, y- 
and z-coordinate separately. Therefore, rotations of simulated bodies cannot be 
modelled with Chain-mail. 
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Schill[10] suggests an enhanced version of Chain-mail (ECM) that processes 
the chains in order of their length violation. This enables ECM to deal with 
inhomogeneous material. 



Our Approaches 

In the following, we suggest three different deformation algorithms that can be 
used separately or combined. 



Quasi-static Euler Modification 



For the simulation of the elastic behaviour of the ILM, we use a modified Mass- 
spring algorithm. 

In a Mass-spring system a high damping can be applied to R^^v) in equa- 
tion 3 to remove energy from the simulated system and increase the stability of 
the simulation. However, a high damping also reduces convergence speed. As- 
suming the simulation is in a steady state after each time-step, the velocity must 
be zero. Setting Ujo_i) = 0 in equation 2 leads to the quasi-static modification 
of the Euler integration: 






mi 



(4) 



Directed Length-Correction 

We suggest a directed version of the already mentioned length correction by 
Provot[6]. The order in which springs are processed depends on the distance to 
the interaction point. Springs near the interaction points are corrected first, far 
springs are corrected later. After altering the position of a node the node is fixed 
for the rest of the correction pass to assure the correction process is directed 
from the interaction point to the edges of the mesh. The algorithm 1 uses a list 
sorted by the distance to the interaction node. 



Algorithm 1 Directed length-correction algorithm. 
1: create sorted spring list 
2: for all strings in list do 
3: correct string length 

4: set string nodes fix for the rest of the pass 

5: end for 



Dragnet Algorithm 

We propose an algorithm called Dragnet to deal with large local displacements. 
The basic idea of Dragnet is pulling on a web of interwoven strings. The point 
where two strings are knotted is called node. Pulling on a node (the interaction 
node) stretches the connected strings. Further dragging of the node leads to a 







